A model of spatiotemporal signal processing by primate cones and horizontal cells

Abstract

A model of spatiotemporal signal processing by the cone–horizontal cell circuit in the primate outer retina is developed and validated using measurements on the H1 horizontal cell from the literature. The model extends an earlier temporal model that mainly addressed the regulation of sensitivity by the cones. Three elements are added to the earlier model to describe the full spatiotemporal processing by horizontal cells. First, the feedback gain from horizontal cells to cones is made adaptive, depending on field size. Second, the spatial filtering by the horizontal dendritic tree is modeled as a two-component spatial filter. Third, an adaptive temporal low-pass filter is added, also depending on field size. The resulting model adequately describes all available measurements on spatiotemporal processing in macaque H1 cells. The adaptive feedback gain is argued to contribute to negative afterimages and chromatic adaptation in human vision.

Introduction

In the primate retina, the transmitter release of cones that drives horizontal and bipolar cells is regulated by feedback from horizontal cells (Wässle, 2004). The H1 horizontal cell receives input from and returns feedback to the medium- and long-wavelength cones exclusively (Dacey, Lee, Stafford, Pokorny, & Smith, 1996). The H1 cell has recently been extensively studied by intracellular recording from isolated macaque retinas, in a preparation with attached pigment epithelium and choroid (Lee, Dacey, Smith, & Pokorny, 2003; Packer & Dacey, 2005; Smith, Pokorny, Lee, & Dacey, 2001). van Hateren (2005) showed that a range of these measurements is well fitted by a model that has two crucial components, namely, phototransduction in the cones and a cone–horizontal cell feedback circuit. Sensitivity regulation and control of integration time in H1 cells are consistent with those expected from the phototransduction cascade, over a range of background intensities spanning three orders of magnitude. The phototransduction part of the model was shown more recently (van Hateren & Lamb, 2006) to be consistent with the photocurrent of human cones as determined from electroretinograms (Friedburg, Allen, Mason, & Lamb, 2004). However, the H1 model was only validated for stimuli of a wide, fixed field size, and it does not address the mechanisms of spatial integration. 

Here, the previous temporal model is extended to include spatial integration. Several additions had to be made to the previous model to be able to describe the measurements adequately. The most important of these is an adaptive feedback gain acting locally at each synaptic triad. A primate cone typically contains a few dozen triads, where each triad consists of two horizontal cell dendrites and one bipolar cell dendrite, positioned in one of the invaginations of the pedicle membrane where transmitter is released (Wässle, 2004). The adaptive feedback gain compensates for deviant transmitter output of the triads by responding with a time constant of at least several seconds. It is argued that this mechanism affects human perception by contributing to negative afterimages and chromatic adaptation. 

Model

The model presented here is developed based on three main requirements, affecting the form of the model and the level of detail chosen. The first and foremost requirement is that it should describe the responses in the outer retina (i.e., cones, horizontal cells, and cone output to the bipolar cells) simply and accurately enough to allow an analysis of the function of the various mechanisms for processing natural stimuli. Examples of how this can be applied can be found in the analysis (van Hateren, 2005, supplementary material) of how the phototransduction cascade contributes to Weber's law (contrast constancy) and how the cones compress the large dynamic range of natural time series of intensities (van Hateren & Snippe, 2006). 

A second requirement of the model is that it should be well posed to act as a preprocessing module for analyzing and modeling downstream neurons. This requires a model that is as simple as possible and that should be in a form that allows highly efficient implementations (making it also feasible to fit the model to measurements). This is accomplished in van Hateren (2005) and here by representing the processing as a simple block scheme, consisting only of first-order low-pass filters and static nonlinearities, and using very fast autoregressive temporal and spatial filters for the computation. 

A final requirement of the model is that it should be as physiologically realistic as possible under the constraint of simplicity; that is, all components should have a clear physiological interpretation. I believe that models developed under such a requirement not only are more likely to prevail in the long run but also stand a better chance to bear fruit as a guide for developing viable full physiological models. The present model is clearly not intended as a full physiological model because many of its components are deliberately abstracted from specific physiological mechanisms. However, I believe that it can serve as a backbone for more elaborated versions once more of the physiological processes in the primate outer retina become firmly established. 

For describing the model, the conventions introduced in van Hateren (2005) will be used. In particular, first-order low-pass filters of unit DC gain are symbolized by a box containing a τ. The gains associated with the physiological substrate of these filters are combined as much as possible with other gains to minimize the number of free parameters. Dimensional conversion constants are omitted to keep the equations as concise as possible. Throughout the article, the term (light) intensity is used when referring to the retinal illuminance expressed in trolands. 

Cones

Figure 1A shows the model for primate cones that has been developed in van Hateren (2005). Light I excites cone pigment R, producing R*. R*, which has a lifetime τ_R, excites a G protein (transducin), which rapidly forms an active complex with phosphodiesterase, yielding E*. E*, which has a lifetime τ_E, hydrolyzes cGMP with a rate β, with β consisting of a dark rate c_β in addition to a factor proportional to E*. An increase in E* thus produces a decrease in cGMP, closing ion channels in the cone outer segment that were held open by cGMP. The current flowing into the outer segment, I_os, subsequently drops. Because part of this current consists of Ca²⁺, the Ca²⁺ concentration in the cell declines, which reduces the Ca²⁺ inhibition of guanylylcyclase to produce cGMP. The original decrease in cGMP is thus counteracted, and these processes form, therefore, a negative feedback loop. 

Figure 1

 

(A) Model of the primate cone. See the text and van Hateren (2005) and van Hateren and Lamb (2006) for details. (B) Horizontal cell feedback model used in the study of van Hateren to describe the temporal responses of macaque H1 cells with stimuli of fixed 5° field size. (C) The extended model adds three new processes to the model of Panel B: The temporal filter τ_p and the feedback gain g_h depend on the strength of the transmitter release I_t, which indirectly depends on the spatial configuration of the stimulus. Furthermore, spatial filtering attributable to the H1 dendritic tree and electrical coupling between H1 cells is added. The inset shows a typical synaptic activation function, Equation 1, with generic parameters. bc = signal toward bipolar cells, hc = measured horizontal cell response.

In van Hateren (2005), it was shown that the nonlinear differential equation that describes the hydrolysis of cGMP can be conceptually and computationally represented by the middle part of the scheme in Figure 1A. The regulation of the gain is dominated by a static nonlinearity 1/β, which is followed by a low-pass filter with time constant τ_X = 1/β. The 1/β nonlinearity and τ_X are directly dependent on the level of E* and, thus, on light intensity, and they are the main factors responsible for the fact that cones become less sensitive and faster when the ambient light intensity becomes higher (Nikonov, Lamb, & Pugh, 2000; van Hateren, 2005). In addition, the calcium feedback has an important role in fine-tuning sensitivity and speed (see the supplementary material in van Hateren, 2005, for an extensive analysis). 

The cone inner segment and pedicle contain voltage- and calcium-sensitive ion channels (Yagi & Macleish, 1994) that shape the properties of the voltage signal arriving at the synapse. In particular, the hyperpolarization-activated current I_h (Hestrin, 1987) produces band-pass filtering (guinea pig rods: Demontis, Longoni, Barcaro, & Cervetto, 1999; modeled cat cones: Hennig, Funke, & Wörgötter, 2002; tiger salamander rods: Mao, MacLeish, & Victor, 2003). Although the influence of these ion channels on macaque cones under physiological conditions is not well known, I assume here (as in van Hateren, 2005) that they are responsible for the low-frequency falloff of the temporal frequency curves (e.g., Figure 5; see further van Hateren, 2005). The analysis in the appendix of Detwiler, Hodgkin, and McNaughton (1980) is adopted to avoid introducing many poorly known parameters, which leads to the rightmost feedback loop in Figure 1A. This loop represents the voltage-sensitive and filtering properties of the inner segment, axon, and cone pedicle (van Hateren, 2005). The resulting membrane voltage arriving at the synapse, V_is, is slightly band-pass filtered by this feedback loop compared with the photocurrent generated by the outer segment (I_os). 

Horizontal cell feedback

Figure 1B shows the model used in van Hateren (2005) to explain responses of the macaque H1 horizontal cell to wide-field stimuli. The horizontal cell feeds back onto the cone pedicle in such a way that its membrane potential, V_h, is subtracted from the membrane potential of the cone, V_is. The result, V_s, effectively drives the synapse through a nonlinear function I_t(V_s), representing the activation function of the synapse. Subsequent low-pass filters τ₁, τ₂, and τ_h then produce V_h. The factor a_I (Figure 1B; see further van Hateren, 2005) is a factor that slightly adjusts two of the time constants and the gain of I_t over the range of background intensities (1–1,000 td) for which the model was validated. For the results presented in this article, with all measurements obtained at high background intensities, a_I was kept fixed at a value of 1. 

The main effect of the subtractive feedback loop of Figure 1B is to suppress low temporal frequencies and enhance high frequencies up to 30–40 Hz. Even higher frequencies are attenuated compared with those in the cone. The overall effect of this is that the output of the cone synapse, feeding into both the horizontal cells and the bipolar cells, becomes band-pass filtered compared with the cones. 

The relatively simple feedback model of Figure 1B works well for wide-field stimuli over a wide range of background intensities and for a variety of temporal stimuli (van Hateren, 2005). However, I found that this model needed three additions to be able to adequately explain the H1 responses when the spatial layout of the stimulus is varied, including narrow-field stimuli with a dark surround. Furthermore, two small changes were made regarding the zero point of the voltages and the equation describing the activation function. I will now first discuss these two changes and then introduce the model additions before discussing them in detail in the next sections. 

A first change is merely administrative. In van Hateren (2005), V_is is defined relative to the voltage that corresponds to I_os = 0. This was convenient for describing the nonlinear dynamics of the feedback loop of the inner segment. However, for describing the feedback loop involving the horizontal cell, in particular the elaborated loop presented here, it is more convenient to define the membrane voltages of cone pedicle and horizontal cell relative to their dark values. Hence, V_iz = V_is − V_is,dark is the voltage of the cone pedicle thus defined. In the dark, V_iz = 0 by definition, and the horizontal cell voltage, V_h, is then zero as well. The latter is guaranteed by choosing the synaptic activation function, I_t(V_s), such that I_t(0) = 0. This choice can be made without loss of generality because it only results in redefining the zero points of the currents and voltages involved. Inspection of the feedback loop in Figure 1C shows that all voltages and currents in the loop must be at their zero points when V_iz = 0. 

The synaptic activation function is chosen here slightly differently from the Boltzmann function used in van Hateren (2005). This is done for mathematical convenience only; the two functions have nearly identical shapes in the region of interest for the stimuli considered here. The activation function is defined as 

$I_{\mathrm{t}}\left(V_{\mathrm{s}}\right)=g_{\mathrm{t}}{V}_{\mathrm{n}}\frac{(s+1)\left(\exp \right(V_{\mathrm{s}}/V_{\mathrm{n}})-1)}{s\exp (V_{\mathrm{s}}/V_{\mathrm{n}})+1},$

(1)

where V_n determines the width of the midsection of the curve, s is a factor that determines the skewness (asymmetry) of the function relative to I_t = 0, and g_t is the dark gain, that is, the slope of the function at V_s = 0. Note that I_t(0) = 0, as required above. The panel in the lower left of Figure 1C shows I_t(V_s) for typical parameter values. V_s will be mostly less than zero, apart from transients, strong surrounds, or strong chromaticity. The activation is therefore mainly situated in the toe of the curve. 

The three processes added to the original model are shown in red in Figure 1C. Firstly, the feedback gain, g _h, is assumed to be adaptive. Secondly, the spatial integration by the horizontal cell is described by two components at different spatial scales (Packer & Dacey, 2002, 2005). Finally, there is an adaptive time constant, τ_p, immediately before the feedback loop. I will discuss below these additions in more detail and justify their inclusion into the model. This is done by changing or omitting, in turn, each of the three additions from the full model of Figure 1C. Thus, the effect of each addition is analyzed, with the other two additions being fully operational. 

Adaptive feedback gain

The symbols in Figure 2A (from Smith et al., 2001) show the sensitivity of a horizontal cell to temporal frequencies produced by modulating a homogeneous stimulus field (mean intensity = 1,000 td, dark surround). The sensitivity for a wide field (10° diameter, open symbols) is two to three times larger than the sensitivity for a narrower field (2° diameter, filled symbols). For comparison, the response of the cone predicted by the model in van Hateren (2005) is shown as the continuous red trace. The shape of this curve deviates strongly from the horizontal cell curve at frequencies just below and above 30–40 Hz. This was modeled in the study of van Hateren as the feedback circuit shown in Figure 1B. The feedback leads to a resonance around 30–40 Hz because the phase delay of the feedback branch is such that those frequencies are approximately 180° out of phase with the forward signal at the point where they are subtracted from that signal. The response is therefore enhanced at those frequencies. At higher frequencies, the feedback gradually ceases to function because the feedback signals are attenuated by the low-pass filters in the loop. As a result, these low-pass filters merely add to the filtering performed by the cone, and therefore, the high-frequency asymptote measured in the horizontal cell is steeper than in the cone. The high-frequency slope is caused by the low-pass filters of the feedback loop (τ₁, τ₂, and τ_h) in addition to the low-pass filters of the cone (van Hateren, 2005). 

Figure 2

 

(A) Black symbols: measurements (from Smith et al., 2001) of an H1 cell's response to temporal sinusoids given on fields with two different diameters. The continuous red line shows the response of the model cone to such modulation. The dashed blue lines show a fit (simultaneously for the amplitude and phase curves at 2°, 5°, and 10° stimulus diameters) according to the model of Figure 1C, with g_h and τ_p fixed. (B) Response of the H1 cell according to the model of Figure 1B, depending on the total gain in the feedback loop. (C) Detail of the circuit involved in adapting the feedback gain g_h. The curve shows Equation 2, with generic parameters. (D) Calculation with the full model of Figure 1C showing the effect on V_d of fixing g_h to different values. The measured V_h (left panel) is rather different for different field sizes because of the cable properties of the H1 cell. The adaptive change in g_h counteracts this, resulting in a V_d that is rather similar at low frequencies (right panel).

Although the model of Figure 1B is adequate for wide-field stimuli, it is confronted with a problem for small-field stimuli. The reduction in amplitude when going from 10° to 2° must be related to the cable properties of the horizontal cell. Horizontal cells form a dense, coupled network of dendrites and, therefore, approximately function as a two-dimensional resistive slab, characterized by a space constant λ (Jack, Noble, & Tsien, 1983). Although the coupling and, therefore, λ depend on the background light level in many vertebrate retinas, this appears not to be the case for macaque H1 cells (Packer & Dacey, 2002, their Figure 11). I will therefore assume that λ is a constant. The horizontal cell network then acts as a fixed spatial low-pass filter on its input (the two-dimensional array of cones). Stimuli of small spatial extent are thus more reduced in amplitude than wide stimuli because of spatial blurring. 

One might think that this reduction in amplitude for small stimuli could fully explain the one observed in Figure 2A when the 10° data are compared with the 2° data. However, this explanation also predicts a strong change in the shape of the sensitivity curve, which is not observed. Figure 2B illustrates this for the model of Figure 1B. It shows the predicted sensitivity curve of the horizontal cell for a range of loop gains, with the loop gain defined as the total gain encountered when going around the entire feedback loop. Part of this gain is determined by the spatial filtering performed by the horizontal cell network. Small-field stimuli therefore produce smaller loop gains than wide-field stimuli. As can be seen in Figure 2B, the sensitivity as observed in the horizontal cell decreases when the loop gain is reduced, which is thus consistent with the differences in sensitivity shown in Figure 2A. However, Figure 2B also shows that as a result of the gain reduction, the resonance is reduced as well. Smith et al. (2001) report that in most cells, the resonance tends to become somewhat smaller when the field diameter is reduced, but the strong reduction seen in Figure 2B is not consistent with the measurements (cf. Figure 2A). Indeed, it can be shown theoretically for the feedback circuit of Figure 1B that a moderate reduction of the amplitude sensitivity at low frequencies necessarily leads to a strong reduction of the resonance, independent of the parameter values used in the model calculation. The dashed blue lines in Figure 2A show the best fit that can be obtained with a feedback model that includes a realistic spatial filter (the two-component spatial receptive field included in Figure 1C) but with a fixed feedback gain g_h. As is clear from the figure, such a model cannot simultaneously fit the measurements at 2° and 10° stimulus diameters. Moreover, the phase characteristic of the 2° model fit deviates strongly from the measured phase in the 20- to 40-Hz region (not shown). 

The only way to obtain a sensitivity reduction at low frequencies while retaining much of the resonance is by decoupling them in the model, that is, by making them (at least partly) depend on different mechanisms. I investigated a range of model configurations to accomplish this, but I found that most configurations can only produce results that are inconsistent with the experimental data. For example, neither changing the shape of the synaptic activation function nor introducing a static nonlinearity for the feedback gain works. The former does not solve the problem illustrated in Figure 2B because it can only change the loop gain and, thus, leaves the low-frequency sensitivity and resonance coupled. The latter is incompatible with the spot-annulus experiment of Lee, Dacey, Smith, and Pokorny (1999; see also below). In this experiment, it was shown that the gain is primarily changed by local changes in illumination and is hardly influenced by modulating the H1 membrane potential V_h by changing the surround illumination. A static nonlinearity for the feedback gain (g_h) driven by V_h can therefore be ruled out as an explanation for the response characteristics discussed above. 

I found only one scheme that could satisfactorily fit the data without becoming unwieldily complex. This scheme is illustrated in Figure 2C. It assumes that g _h, the feedback gain exerted by each horizontal cell dendrite onto the presynaptic elements in its own triad, is not fixed but adaptive (see the Discussion section for possible physiological substrates). The measured membrane potential of the horizontal cell, V _h, is then different from the effective signal produced by the dendrite, V _d, which is the actual signal subtracted from the potential at the cone pedicle, V _p. The result, V _s, drives the transmitter release of the triad. Depending on the physiological substrate of the feedback mechanism (which is still uncertain, see the Discussion section), V _d may be interpreted as, for example, the extracellular voltage surrounding the dendrite or the local H ⁺ concentration. 

Note that this model assumes that the feedback and the regulation of its gain are strictly local to the triad: They only depend on the voltage in the horizontal cell dendrites in the triad and on the pedicle voltage. The feedback is therefore not directly coupled to the other triads of the same pedicle. This lack of feedback coupling will produce the observed lack of (spectral) cross talk between H1 and H2 horizontal cells (Dacey, 1999) if the horizontal cell dendrites within a particular triad are preferentially of similar type (both are H1 or both are H2). 

I found that the data of Smith et al. (2001) could be well fitted by assuming that g_h is driven by the output I_t of the synaptic activation function of the triad, via a low-pass filter τ_itd and a nonlinearity NL_itd. The time constant τ_itd must be longer than a few seconds because, otherwise, its effect would have been detected (as gain changes dependent on V_h) in the spot-annulus experiment of Lee et al. (1999). These experiments, with an annulus modulating V_h (and thus, indirectly, I_t) at 0.6 Hz, showed no corresponding gain changes for stimuli presented in the central spot, and therefore, τ_itd must be significantly longer than 1/0.6 s. Simulations confirm that a τ_itd shorter than a few seconds does not produce adequate fits to these measurements (see further the Spot-annulus experiment section). 

There is no obvious upper limit for τ _itd because the experiments of Smith et al. (2001) and Packer and Dacey (2002, 2005) were all done with a static field diameter and, therefore, provide no specific information on τ_itd. I will assume τ_itd = 10 s below, a value that is sufficiently long to obtain steady-state conditions for the measurements considered here but is otherwise rather arbitrary. 

The nonlinearity NL _itd is also poorly constrained because only a few field diameters were used in the experiments (typically 2°, 5°, and 10°). I found that a simple nonlinearity of the form  

$g_{\mathrm{h}}\left(I_{\mathrm{t}}\right)=\frac{1}{1+\exp \left[c_{\mathrm{h}}\right(I_{\mathrm{t}}-I_{\mathrm{h}}\left)\right]},$

(2)

where c _h and I _h are constants, was able to adequately describe the available measurements with different field sizes. The argument of g _h is here formally I _t, as it would be in the steady state (when the τ _itd filter has unit gain), but more generally, its argument would be the result of I _t low-pass filtered by τ _itd. An example of g _h with typical parameter values is shown in Figure 2C. 

The effect of a variable g _h is illustrated in Figure 2D, which shows calculations with the full model of Figure 1C, but with g _h clamped to particular values. Whereas the V _h for 10° and 2° field diameters are quite different at low frequencies because of the cable properties of the horizontal cell, this is partly compensated in the feedback signal V _d by an increased g _h for the smaller field diameter. The result is that the values of V _d for 10° and 2° field diameters are not very different at low frequencies, and the loop gain for 2° is less reduced compared with the gain for 10° than would have been the case without the change in g _h. Consequently, V _h for 2° still shows a considerable enhancement of frequencies up to 30–40 Hz compared with the cone response, in accordance with the measurements on H1 cells. 

A two-component spatial receptive field

Recently, it was reported (Packer & Dacey, 2002, 2005) that macaque H1 cells can be well described by a spatial receptive field consisting of two components acting on different spatial scales. The component with the shorter scale, λ_S, is presumably mainly determined by the H1 cell's own dendritic tree, whereas the component with the longer scale, λ_L, presumably reflects the electrical coupling between H1 cells (Packer & Dacey, 2005). I found that it is indeed necessary to include such a two-component spatial receptive field into the model to explain the measured spatial frequency response. 

Figure 3A illustrates the problem when the spatial integration of the H1 cell is represented by a single space constant λ (according to the left diagram in Figure 3B). The best fit with all other components of the model of Figure 1C included shows a clear discrepancy with the measurements (red trace: fit; symbols: measurements from Packer & Dacey, 2002). For the sake of an efficient computation (see the Model implementation section), the spatial filtering was implemented as a filter separable from the temporal filtering of the H1 cell. For the two-dimensional cable equation, this is quite accurate on a time scale slower than the time constant of the H1 cell but is only an approximation on faster time scales. Provisional calculations using the two-dimensional cable equation indicate that the lack of spatiotemporal inseparability does not affect the conclusions presented in this article. 

Figure 3

 

(A) Measurements of the spatial frequency response (black symbols: data from the upper trace of Figure 4 in Packer & Dacey, 2002) and a fit (red line) using a single length constant. The inset shows the same data in a double logarithmic plot. (B) The single length constant (left diagram) is replaced by a weighted superposition of filters with short and long length constants (right diagram). (C) The fit resulting from the two-component spatial filter. See text for further details. Spatial frequencies are given in cycles per micrometer; this may be converted to cycles per degree by assuming that, in the macaque eye, a 200-μm retinal distance corresponds approximately to 1° of visual angle (Smith et al., 2001).

The two-dimensional spatial filter is implemented here as an autoregressive filter (see the Model implementation section) acting on a hexagonal sampling grid. The spatial point spread function of this filter, h( r), approximates an exponential; that is,  

$h\left(r\right)\approx c_{\lambda }\exp (-r/\lambda ),$

(3)

where r is the distance from the origin, λ is the space constant, and c _λ is a normalization constant such that the filter has unit response to a homogeneous field. It should be noted that the two-dimensional cable equation has a modified Bessel function K ₀ as its (steady-state) spatial point spread function (Jack et al., 1983) and that an exponential only approximates this at large distances. Whereas the K₀ function gives 1/(1 + (2πf_sλ)²) as the corresponding spatial modulation transfer function, with f_s spatial frequency, an exponential has 1/(1 + (2πf_sλ)³) as its spatial modulation transfer function. The analysis performed in Packer and Dacey (2002) on the spatial frequency curves used the latter. However, I found that performing the same analysis on the same data gives rather similar results for either choice of filter: In both cases, a superposition of two filters is needed to explain the measurements. Here, I will use filters based on Equation 3. 

The diagram at the right in Figure 3B shows how the two spatial scales were implemented as a superposition, giving a point spread function  

$h\left(r\right)=w_{\mathrm{S}}{c}_{{\lambda }_{\mathrm{S}}}\exp (-r/{\lambda }_{\mathrm{S}})+(1-w_{\mathrm{S}})c_{{\lambda }_{\mathrm{L}}}\exp (-r/{\lambda }_{\mathrm{L}}),$

(4)

where λ _S ( λ _L) is the short (long) space constant and w _S (with 0 ≤ w _S ≤ 1) is the relative weighting of the λ _S component. Note that the normalization of the two contributing filters ensures that h( r) has a unit response to a homogeneous field. 

Figure 3C shows the fit with this model to the data (red trace and black symbols). Because of the feedback and because of the nonlinearity in the circuit, linear superposition does not hold. It is nevertheless instructive to investigate the approximate contributions of the short- and long-scale components to the spatial frequency response. To avoid shifting the DC level (for f _s = 0) and thereby strongly changing the working point of I _t( V _s), this was done in the following way. First, the calculation was performed assuming only the filter with λ _S (with unit weight, as in Equation 3), while keeping all other parameters fixed to their fitted values. The resulting response was multiplied by w _S and shown as the blue lines marked λ _S in Figure 3C. Second, the calculation was performed assuming only the filter with λ _L, and the result was multiplied by 1 − w _S (blue lines marked λ _L). Although it is clear that there is no linear superposition in the midrange of frequencies, it can be seen that, as expected, λ _L is mainly responsible for the peak of the spatial frequency response at low frequencies and that λ _S is mainly responsible for the tail at high frequencies. 

Adaptive temporal filtering

The final addition to the original model is relatively minor compared with the adaptive feedback gain and the spatial filtering discussed above. Figure 4A shows the high-frequency sections of the sensitivity curves of three H1 cells (a, b, and c; from Smith et al., 2001), measured at field sizes 10° and 2° (black open symbols). The red dots show the fits made with the complete model of Figure 1C but without the branch marked “adaptive temporal filtering.” Although the fits are fairly good, close inspection shows a slight but systematic deviation. The resonance at 10° field size is slightly too strong, and the resonance at 2° field size is slightly too weak. Related to this, the frequency falloff for the 2° field size is too shallow. Corresponding deviations occur in the phase characteristics of the cells (not shown), making the actual fits worse than would appear from Figure 4A. Because these deviations are present in all three cells, they indicate a missing mechanism in the model. The too shallow frequency falloff at 2° field size suggests that there is stronger low-pass filtering at 2° than at 10°. I found that making one or more of the low-pass filters inside the feedback loop adaptive, depending on field size, did not solve the problem. The only simple and effective way to improve the fit was by introducing an adaptive low-pass filter, τ_p, immediately before the feedback loop (see the fits in Figure 4B). The change in τ_p is assumed to be driven by the output I_t of the synaptic activation function, via a slow low-pass filter, τ_itp, and a nonlinearity, NL_itp (Figure 4C). For the same reasons as discussed above for τ_itd, τ_itp must be slower than a few seconds but has no upper bound based on the measurements considered here. I will assume τ_itp = 10 s below. The measurements can be well fitted by assuming a nonlinearity, NL_itp, of the same form as NL_itd: 

${\tau }_{\mathrm{p}}\left(I_{\mathrm{t}}\right)=\frac{{\tau }_{\mathrm{p},\max }}{1+\exp \left[c_{\mathrm{p}}\right(I_{\mathrm{t}}-I_{\mathrm{p}}\left)\right]},$

(5)

where c_p, I_p, and τ_p,max are constants. An example of τ_p with typical parameter values is shown in Figure 4C. For a field size of 2°, V_d is smaller (closer to zero, thus less close to V_iz and V_p) than for a field size of 10° (which is less reduced by the cable properties of the horizontal cell). Therefore, V_s is more negative (closer to V_iz and V_p) for 2° than for 10°; hence, I_t is more negative, and therefore, τ_p is larger for 2° than for 10°. 

Figure 4

 

(A) Measurements from three H1 cells (a, b, and c; from Smith et al., 2001) at 2° and 10° field sizes (black symbols, connected by straight lines for the sake of clarity) and fits with a model without an adaptive low-pass filter τ_p (red dots, connected by straight lines). (B) Same as in Panel A, with fits using the model with adaptive τ_p. (C) Detail of the circuit involved in an adaptive τ_p. The curve shows Equation 5, with generic parameters.

Model implementation

Temporal filtering

All temporal filters in the model are first-order low-pass filters (see Figures 1A and 1C) and were implemented as autoregressive moving-average (ARMA) filters (Brown, 2000; van Hateren, 2005). The output y(n) to an input x(n) is then given by 

$y\left(n\right)=f_{1}y(n-1)+f_{2}x(n-1)+f_{3}x\left(n\right),$

(6)

with 

$\begin{array}{l}{f}_1=\exp (-1/{\mathrm{\tau \prime }}^{\prime })\\ f_2={\mathrm{\tau \prime }}^{\prime }-(1+{\mathrm{\tau \prime }}^{\prime })\exp (-1/{\mathrm{\tau \prime }}^{\prime })\\ f_3=1-{\mathrm{\tau \prime }}^{\prime }+{\mathrm{\tau \prime }}^{\prime }\exp (-1/{\mathrm{\tau \prime }}^{\prime })\\ {\mathrm{\tau \prime }}^{\prime }=\tau /\Delta t\end{array}$

(7)

where τ is the time constant of the low-pass filter and Δt is the time step. For the calculation, Δt = 0.1 ms turned out to be sufficiently short for accurate results. This was verified, firstly, by varying Δt, and, secondly, by comparing the results with those obtained by numerical integration in Matlab. An ARMA implementation (in Fortran) is about two orders of magnitude faster than solving the system of differential equations by conventional numerical integration and, therefore, crucial for obtaining fits within a reasonable amount of time (see the Computation section). 

Spatial filtering

An implementation of the spatial filter that is equally efficient as that of the temporal filter was obtained in the following way. A filter with a one-sided exponential pulse spread function, converting a spatial input p( n) into a spatial output q( n), is given by  

$q\left(n\right)=g_{1}q(n-1)+g_{2}p\left(n\right)$

(8)

with n increasing, and  

$\begin{array}{l}{g}_1=\exp (-1/{\mathrm{\lambda \prime }}^{\prime })\\ g_2=1-\exp (-1/{\mathrm{\lambda \prime }}^{\prime })\\ {\mathrm{\lambda \prime }}^{\prime }=\lambda /\Delta s\end{array}$

(9)

where λ is the space constant of the filter and Δ s is the step size. For most calculations, Δ s = 0.3° (with field size 10°) turned out to be sufficiently small for accurate results; for Figure 6, it was necessary to use Δ s = 0.1° (with field size 15°), and for Figure 8, Δ s = 0.05° (with field size 5°) was needed. The latter values of Δ s are close to the cone spacing expected at the approximate eccentricity of the measured cells. Accuracy was checked by varying Δ s and field size. 

Note that the spatial filter of Equations 8 and 9 is slightly different from the temporal filter defined by Equations 6 and 7. In particular, f ₂ and f ₃ are combined into g ₂, and only the present input sample, p( n), is used and not the previous input sample, p( n − 1). The reason is that the temporal filter of Equations 6 and 7 is carefully crafted (Brown, 2000) to give the correct phase characteristics of the (causal) low-pass filter, given a finite time step Δt. Causality is of no concern when dealing with spatial filters; hence, Equations 8 and 9 provide the optimal filter, with an impulse response that is simply a sampled exponential irrespective of the step size Δs. 

Applying subsequently the same filter into the reverse spatial direction, that is, q( n) ≔ g ₁ q( n + 1) + g ₂ q( n) for decreasing n, turns out to produce a perfect two-sided exponential pulse spread function. In the equation above, ≔ denotes assignment, to indicate that the q( n) to the left is different from the q( n) to the right. For finite arrays [1: N], boundary effects have to be taken into account. It can be shown that by choosing the first output element for the filtering into the forward direction as q(1) = g ₁ b + g ₂ p(1), and the first output element for the filtering into the reverse direction as q( N) ≔ ( q( N) + g ₁ b)/(1 + g ₁), the filter output is identical to that of an infinite input array with the values beyond the boundary given by b (thus, b is the intensity of the surround if the input is an array of intensities). 

The two-sided exponential filter described above was extended to two dimensions by consecutively applying it along the three main axes of a hexagonal (i.e., triangular) grid of sampling points (i.e., cones). Because the three axes are not orthogonal, the filtering operations become coupled. However, the resulting point spread function is still exponential to good approximation. The point spread function is approximately circularly symmetrical because the two-sided exponential filter applied to each axis is linear, and the order of processing, therefore, does not matter (hence, all axes are necessarily equal). The coupling mentioned above produces boundary effects that slightly disturb this equality, but I found that their effect on the present calculations is negligible. 

The full spatial filter was obtained by separately filtering with λ _S and with λ _L and by combining the results according to the weighting of Equation 4. 

Computation

The model was implemented in Fortran and run in parallel on an 8-processor (3 GHz Xeon) Beowulf cluster. Fitting was performed using a simplex algorithm (Press, Teukolsky, Vetterling, & Flannery, 1992) for minimizing the RMS deviation between model responses and measurements. Despite the parallel processing and the highly optimized filter algorithms, simultaneously fitting the responses to many stimuli (e.g., all frequencies and field sizes in Figure 5A or all pulses in Figure 6) typically took several hours. A single model calculation for a 1-s stimulus with Δs = 0.3° and 10° field size (i.e., a circular hexagonal grid with about 1,000 cones) takes 8 s on a single processor. 

Figure 5

 

(A) Temporal sensitivity of an H1 cell for stimulus field sizes of 2°, 5°, and 10°. Black symbols: data from Figure 3A of Smith et al. (2001). The red lines show a fit with parameter values τ_R = 3.7, τ_E = 12, k_β = 4.9 × 10⁻⁵, a_C = 0.30, a_is = 4.8 × 10⁻², t_delay = 2.7, g_t = 21, τ₁ and τ₂ = 6.3, τ_h = 7.0, λ_S = 4.7, λ_L = 450, w_S = 0.22, τ_p,max = 18, c_p = 0.34, I_p = −21, c_h = 0.41, and I_h = −17. Units are given in Table 1, and parameters that were not mentioned are used at their generic values (Table 1). The dashed blue lines show the result using generic values for all parameters. (B) Same as in Panel A (data from Figure 3B of Smith et al., 2001). Parameter values: τ_R = 1.4, τ_E = 19, k_β = 4.7 × 10⁻⁵, a_C = 0.32, a_is = 7.7 × 10⁻², t_delay = 2.5, g_t = 17, τ₁ and τ₂ = 6.9, τ_h = 4.4, λ_S = 22, λ_L = 380, w_S = 0.20, τ_p,max = 32, c_p = 0.46, I_p = −22, c_h = 0.44, and I_h = −18. (C) Same as in Panel A (data from Figure 3C of Smith et al., 2001). Parameter values: τ_R = 3.8, τ_E = 4.9, k_β = 6.2 × 10⁻⁵, a_C = 0.20, a_is = 3.5 × 10⁻², t_delay = 2.8, g_t = 13, τ₁ and τ₂ = 4.1, τ_h = 9.6, λ_S = 22, λ_L = 260, w_S = 0.27, τ_p,max = 20, c_p = 0.22, I_p = −26, c_h = 0.12, and I_h = −21.

Figure 6

 

Responses to 10-ms flashes (denoted by a black bar on the abscissa) given on stimulus fields of 15° (A), 4° (B), 2° (C), and 1° (D), all at 316 td. Flashes have the same extent as the fields. The numbers 1, 2, 4, and 9 at the traces give the amplitude of the flashes in terms of the field intensity (i.e., “2” denotes a 732-td flash given on top of the 316-td field). Responses to Contrasts 2, 4, and 9 have been shifted downward from 0 for the sake of clarity. Black symbols: data from Figure 7 of Smith et al. (2001). Red lines: simultaneous model fit to all responses, with the following parameter values: τ_R = 2.9, τ_E = 5.9, k_β = 3.0 × 10⁻⁵, a_C = 0.13, a_is = 0.22, t_delay = 1.3, g_t = 19, V_n = 8.1, s = 0.32, τ₁ and τ₂ = 6.6, τ_h = 6.0, λ_S = 19, λ_L = 210, w_S = 1.7 × 10⁻³, τ_p,max = 34, c_p = 0.34, I_p = −18, c_h = 0.10, and I_h = −27. Blue lines: response of the generic model.

Comparison with measurements

Fits of the model of Figure 1C to a range of different stimuli will be shown below. Depending on the particular stimulus, the set of parameters to be varied in each fit differs because there are usually some parameters that have no or little influence on the precise response of the model to that particular stimulus. Such parameters were therefore fixed to default (“generic”) values when fitting to that stimulus. Details on the fitted parameter values are given in the figure captions, and the generic parameter values are listed in Table 1. 

Table 1

 

Parameters used in the model (see Figure 1 and Equations 1, 2, 3, 4 and 5). Generic values are used for calculating the generic curves in the figures (dashed blue lines). The range shows the minimum and maximum values obtained from all fits presented here. Notes: the smallest of the values at τ _R and τ _E is arbitrarily assigned to τ _R; for τ _m, the value reported by Hornstein et al. (2004) was used; for simplicity, the time constants τ₁ and τ₂ were taken as equal; t_delay is an overall delay of the response, partly attributable to various diffusional delays distributed over the various model stages. au = arbitrary unit.

Symbol	Description	Units	Generic value	Range
τ R	Time constant of R* inactivation	ms	3.4	1.4 to 3.8
τ E	Time constant of E* inactivation	ms	8.7	4.5 to 19
c β	Rate constant of cGMP hydrolysis in darkness	(ms) −1	2.8×10 −3	Fixed
k β	Rate constant of cGMP hydrolysis	(ms) −1/td	1.0×10 −4	4.7×10 −5 to 3×10 −4
n X	Apparent Hill coefficient of CNG activation	–	1	Fixed
τ C	Time constant of Ca 2+ extrusion	ms	3	Fixed
a C	Scaling constant of GC activation	au	0.2	0.09 to 0.32
n C	Hill coefficient of GC activation	–	4	Fixed
τ m	Capacitive membrane time constant	ms	2.3	Fixed
γ	Parameter of membrane nonlinearity	–	0.7	Fixed
a is	Scaling constant of membrane nonlinearity	au	0.1	3.5×10 −2 to 3.7×10 −1
τ is	Time constant of membrane nonlinearity	ms	90	Fixed
t delay	Delay time	ms	3	1.3 to 3.5
g t	Parameter of transmitter activation curve	au	15	8 to 21
V n	Parameter of transmitter activation curve	mV	5	3.3 to 8.1
s	Parameter of transmitter activation curve	–	0.4	0.3 to 0.4
τ 1, τ 2	Time constants of cone–horizontal cell loop	ms	5.7	4 to 8
τ h	Time constant of cone–horizontal cell loop	ms	7.0	4 to 20
λ S	Short length constant of spatial filter	μm	20	5 to 30
λ L	Long length constant of spatial filter	μm	300	100 to 440
w S	Relative weight of short spatial filter	–	0.15	0 to 0.27
τ p,max	Parameter of adaptive τ p loop	ms	25	18 to 34
c p	Parameter of adaptive τ p loop	au	0.25	0.22 to 0.46
I p	Parameter of adaptive τ p loop	au	−25	−18 to −26
τ itp	Parameter of adaptive τ p loop	s	10	Fixed
c h	Parameter of adaptive g h loop	au	0.25	0.10 to 0.44
I h	Parameter of adaptive g h loop	au	−20	−17 to −27
τ itd	Parameter of adaptive g h loop	s	10	Fixed

Because the model is complex and contains many parameters, there is the potential threat of overfitting; that is, details of the responses may be specifically fitted by unrealistic values (outliers) of some of the parameters. To show that this is not a significant problem for the model presented here, I determined a set of generic parameter values. These values were chosen to lie well within the range of parameter values obtained from the fits and to produce typical H1 responses. These generic responses are, in most figures, shown as dashed blue lines, in addition to the data (black symbols in the figures) and the fits (solid red lines). Although the generic model predictions differ quantitatively from the data and fitted curves shown, they possess all the qualitative characteristics of H1 responses. The fits to individual cells are obtained for parameter values that differ only modestly from the generic values ( Table 1). This would not be the case if the fits were the result of overfitting with an unrealistic model. Instead, it appears that the model contains the minimum of mechanisms needed for describing signal processing by H1 cells. 

A further indication that the model and its parameters are sufficiently constrained by the data is the fact that values obtained here for parameters of the cone part of the model are very similar to those obtained from other measurements (macaque H1 cell: van Hateren, 2005; human ERG: van Hateren & Lamb, 2006). However, a full evaluation of how well the model generalizes across new types of stimuli can only be obtained when it is confronted with new measurements, either from the H1 cell or from downstream neurons driven by the output of the cone–horizontal cell circuit. The latter will be considered in a forthcoming publication. 

Temporal frequency curves

Figure 5 shows the sensitivity and phase characteristics of three H1 cells (A, B, and C; measurements from Smith et al., 2001). The stimulus was a sinusoidally modulated homogeneous field of 10°, 5°, or 2° diameter, with an average intensity of 1,000 td and a dark surround. The measured sensitivity was determined from the first harmonic of the response to low-contrast (up to 25%) modulations (Smith et al., 2001), that is, in a region where the contrast-response characteristic is still linear to good approximation. The computed response (red lines) was similarly obtained from the model response to 25% contrast sinusoidal modulations. 

Both measured and computed responses are smaller for small-diameter stimulus fields than for large-diameter stimulus fields. This difference can be attributed to the spatial modulation transfer function of the H1 cell dendritic network. In addition, the resonance at 30–40 Hz is slightly reduced at small-field diameters (particularly clear in Cells B and C). The resonance is visible either as a peak (Cell A) or as a sharp corner at the transition from the almost-flat low-frequency part of the curves to the steep high-frequency falloff (Cells B and C). The resonance is also reflected in a fast phase change around 30–40 Hz (lower panels), again in both the measured and model responses. The generic curves (dashed blue lines) follow the same general characteristics. 

Pulse responses

The responses to short pulses of various amplitudes superimposed on the full extent of fields of various sizes are shown in Figure 6 (measurements from Smith et al., 2001). The 10-ms pulses (indicated by horizontal bars on the abscissa) start at Time 0 and have Amplitudes 1, 2, 4, and 9 relative to the background. The responses to Contrasts 2, 4, and 9 have been shifted downward along the ordinate in Figure 6 for the sake of clarity. The measured and fitted model responses share the same basic characteristics: The responses become considerably smaller for smaller field diameters (note the difference in scales), resonance becomes weaker for smaller field diameters, and the nonlinear dependence on contrast is the same in measurements and model response (note that the response to Contrast 9 is considerably less than nine times the response to Contrast 1). However, there are also some discrepancies. One is particularly apparent at Contrast 9 for the 15°, 4°, and 2° field sizes, where the rising flank of the model responses is steeper than that of the measured responses. For small contrasts, there is no such difference. The impulse response of a linear system should display increasing slopes of the rising flank for increasing stimulus contrasts. However, close inspection of the measured curves in Figure 6A suggests that the slope of the rising flank does not increase further after reaching a maximum slope; that is, the slope appears to saturate. This is a type of nonlinearity that is well known in electronics, where it is known as the slew rate limit of an amplifier. It appears to correlate here mainly with the stimulus contrast and not with the amplitude of the voltage swing of the H1 cell (cf. 2° Contrast 9 with 15° Contrast 1, with similar voltage swings but rather different slew rates). It is therefore most likely attributable to individual cones or synapses. Because no slew rate limit is present in ERG-derived photocurrent responses of human cones to high-contrast flashes (Friedburg et al., 2004; van Hateren & Lamb, 2006), it probably originates from processes occurring in the cone pedicle or synaptic triad, such as the mechanism of vesicle release. I decided not to include a specific slew rate limit into the present model firstly because the discrepancy only occurs for rather high contrasts with a fairly low rate of occurrence in natural stimuli and, secondly, because it may be more adequately modeled once more details are known on the dynamical interactions of the processes in the cone pedicle and synapse. 

A second discrepancy can be seen in Figure 6D, where the fitted responses have a somewhat longer latency than the measured responses. This may arise from the choice, made for reasons of computational efficiency, to model the spatial processing of the H1 cell as separable from the temporal processing. In reality, small stimuli are expected to be processed with a shorter latency by the H1 cell because the effective time constant ( RC time) of the cell is expected to be shorter for small fields than for large fields. This is expected because current can also flow laterally in the case of small fields (making R and, thus, the time constant effectively smaller) but less so in the case of large fields when the cell is more isopotential. It is possible that including this into the model would resolve the discrepancy in Figure 6D. 

The curves calculated for the generic model parameters follow the general trend. It should be understood that the precise shape of the responses depends strongly on the spatial characteristics of the particular H1 cell ( λ _S, λ _L, and w _S), and these can be quite different for different H1 cells even at a similar eccentricity (Packer & Dacey, 2002, 2005). 

Spot-annulus experiment

In Lee et al. (1999), a spot-annulus experiment has been presented to show that a considerable part of the gain control as observed in the H1 cell is attributable to a spatially localized process. Figure 7 shows that the present model is well capable of replicating the experimental results. The upper panel shows the experiment where a low-frequency, high-contrast vehicle stimulus (mean intensity, 1,000 td; frequency, 0.61 Hz; contrast, 90%) was given to an annulus (inner diameter, 1.5°; outer diameter, 10°), resulting in a slow modulation of the H1 membrane potential. Simultaneously, a high-frequency, low-contrast test stimulus was given to the central spot (filling the hole in the annulus; diameter, 1.5°; mean intensity, 1,000 td; frequency, 19.52 Hz; contrast, 10%), resulting in a fast response modulation superimposed on the vehicle response (Figure 7A). 

Figure 7

 

(A) Black symbols: measurements from Figure 4 of Lee et al. (1999), where an annulus was modulated at a vehicle frequency of 0.61 Hz (slow wave) and the central spot was modulated at a test frequency of 19.52 Hz (fast ripple). See the text for further explanation. Red line: fitted parameters: τ_E = 4.5, k_β = 5.0 × 10⁻⁵, a_C = 0.18, a_is = 0.11, t_delay = 3.5, λ_S = 30, λ_L = 420, and w_S = 0.15; the adaptive feedback loops were kept fixed, with τ_p = 3 and g_h = 0.5. (B) Same as in Panel A, with spatially superimposed vehicle and test waves. Parameter values are the same as in Panel A. Blue lines: generic model.

If the gain at the central spot is not modulated by the vehicle stimulus given on the surround, the response to the test modulation should have the same amplitude irrespective of its position in the vehicle wave. This is indeed approximately the case in Figure 7A, in both the measured and fitted responses. However, close inspection reveals a slight change in amplitude, with the test response slightly smaller when the vehicle response is positive (i.e., when the annulus is darkest) than when it is negative. This effect is also reported by Lee et al. (1999). The origin in the model is the nonlinear synaptic activation function (Figure 1C): When the H1 response V_h (and consequently V_d) becomes more positive, V_s becomes more negative, thus shifting the working point to shallower parts of the synaptic activation function. As a result, the modulation of V_s attributable to the test stimulus produces smaller modulations in I_t than when V_h is negative. The effect is rather modest in the fitted curve in Figure 7A but larger in the generic curve (dashed blue line). For different H1 cells, the effect varies from small to considerable (B. B. Lee, personal communication), which may correspond in the present model to slight variations in the shape of the synaptic activation function and in the working point. Kraaij, Spekreijse, and Kamermans (2000) have investigated the influence of the nonlinear activation function on the feedback in fish horizontal cells, and the consequences for color vision are discussed in Kamermans, Kraaij, and Spekreijse (1998). 

In Figure 7B, the vehicle modulation was shifted from the annulus to the central spot and, thus, given spatially superimposed on the test modulation. The contrast of the vehicle was reduced to 60% to obtain approximately the same voltage swing in the H1 cell as before. Now, there is a strong modulation of the test response, depending on the intensity of the vehicle. When the vehicle is dark (positive H1 potential), the test response is much larger than when the vehicle is bright. This implies that the test response is larger when the test contrast (test intensity modulation divided by momentary vehicle intensity) is large. This gain control was shown in van Hateren (2005) to be a property of the cone part of the model, mainly attributable to the 1/β nonlinearity in combination with the calcium feedback loop (Figure 1A). 

In principle, the modulation of the feedback gain ( Figures 1C and 2C) might influence the present model fits if the time constant τ _itd is not as long as assumed (10 s). Therefore, I checked if a short τ _itd would be consistent with the measurements in Figure 7. I found that the fits deteriorated when τ _itd was made small enough to influence the dynamics on the time scale of the experiment (∼1.6 s). This supports the notion that the alleged adjustment of the feedback gain is indeed a rather slow process that is not easily detected with stimuli confined to relatively short time scales. 

Spatial frequency curves

Figure 8 shows the spatial frequency response of four different H1 cells (data from Packer & Dacey, 2002). The cell shown in Figure 8A is the same as that in Figure 3C and is reproduced here to facilitate comparison with the other cells and with the generic response (dashed blue lines). As reported by Packer and Dacey (2002, 2005), the spatial frequency response of H1 cells varies considerably from cell to cell, not only as a function of eccentricity but also at a fixed eccentricity. On a linear scale, there is usually a narrow peak at low spatial frequencies and a long tail at high spatial frequencies. Often, there is a notch at medium frequencies (Packer & Dacey, 2002, 2005). The narrow low-frequency peak corresponds in the space domain to a long tail in the point spread function of the H1 cell, presumably attributable to the extensive electrical coupling between H1 cells (Packer & Dacey, 2002, 2005). The long high-frequency tail corresponds in the space domain to a fairly sharp peak in the point spread function of the H1 cell (Packer & Dacey, 2002), presumably attributable to the cell's own dendritic tree. Figure 8 shows that the two-component receptive field used in the model provides a good description of the spatial frequency curves of all cells shown. 

Figure 8

 

Spatial frequency response: data from the upper four curves of Figure 4 in Packer and Dacey (2002). Insets show the same data in double logarithmic plots. Parameter values: (A) a_is = 0.14, λ_S = 18, λ_L = 96, and w_S = 0.13; (B) a_is = 0.37, λ_S = 22, λ_L = 140, and w_S = 0.11; (C) a_is = 4.2 × 10⁻², λ_S = 12, λ_L = 160, and w_S = 4.8 × 10⁻²; (D) a_is = 0.29, λ_S = 24, λ_L = 230, and w_S = 3.3 × 10⁻². The adaptive feedback loops were kept fixed at Panels A–D, with τ_p = 3 and g_h = 0.5. Blue lines: generic model. The thin black line in the inset of Panel A shows the effect of a 10-fold reduction of the feedback gain.

The cone–horizontal cell feedback loop reduces low spatial frequencies more than high spatial frequencies. To illustrate the effect of feedback, the thin black line in the inset of Figure 8A shows the spatial frequency response with the feedback gain reduced 10-fold. High spatial frequencies are little affected by the strength of the feedback, whereas low spatial frequencies are reduced when the feedback becomes stronger. Because this also affects the inferred space constant of the horizontal cell network, this effect needs to be taken into account when a detailed comparison of physiological and dendritic receptive fields is made. 

Effects of carbenoxolone and cobalt

Applying carbenoxolone to the retina interferes with the feedback from horizontal cells to cones, presumably by blocking hemichannels in the horizontal cell dendrites (Kamermans et al., 2001). This effect occurs not only in goldfish (Kamermans et al., 2001) and turtle (Pottek et al., 2003) but also in macaque H1 cells (Packer & Dacey, 2005) and in feedback signals observed in macaque cones (Verweij, Hornstein, & Schnapf, 2003). Similar effects are reported for Co²⁺ (Fahrenfort, Sjoerdsma, Ripps, & Kamermans, 2004; Packer & Dacey, 2005). The main effect is that the slow sagging in the response to an on-step is strongly reduced and that the response amplitude to such steps gradually decreases in the course of the application of the drug and eventually vanishes altogether. 

The large changes caused by carbenoxolone are likely to bring the system far beyond the normal working range assumed for the present model. For example, the ionic driving forces will change considerably, which is at a level of physiological realism beyond the scope of the model. Nevertheless, the model is, in principle, capable of mimicking the effects of carbenoxolone (and similarly of Co ²⁺), as shown in Figure 9. The black traces at Δ V _s = 0 mV and at Δ V _s = −6.5 mV are measured responses (Packer & Dacey, 2005) to a 100% modulated square wave (2.44 Hz, 10° field diameter, 1,000 td mean intensity). The curve at ΔV_s = 0 mV is the response as measured before application of carbenoxolone, and the response at ΔV_s = −6.5 was recorded a few minutes after the start of carbenoxolone application, when the response dynamics were clearly affected, but before the response had vanished (Packer & Dacey, 2005). The red traces are model calculations of the H1 response, V_h, made on the assumption that the main effect of carbenoxolone (or Co²⁺) is to gradually shift V_s (see Figure 1C) to more negative values (denoted by the numbers, ΔV_s in millivolts, below the curves). Two equivalent ways to put this are, first, that carbenoxolone gradually shifts the synaptic activation function to more positive values or, second, that carbenoxolone gradually shifts the feedback signal V_d to more positive values. The model predictions based on this assumption are shown by the red traces, which fit well with the measured traces at ΔV_s = 0 mV and at ΔV_s = −6.5 mV. The slow sagging after an on-step gradually disappears, whereas the slow sagging after an off-step survives longer. Moreover, the response eventually vanishes altogether. 

Figure 9

 

Measured (black lines: first two traces from Cell 1 of Figure 7 of Packer & Dacey, 2005) and computed (red lines) responses before and during application of carbenoxolone. See text for further explanation. Parameter values: k_β = 1.6 × 10⁻⁴, a_C = 0.09, a_is = 0.24, g_t = 8, V_n = 3.3, s = 0.3, τ₁ and τ₂ = 8, τ_h = 20, λ_L = 240, and w_S = 0.25. For simplicity, the adaptive feedback loops were kept fixed, with τ_p = 15 and g_h = 0.5.

The way this works in the model is illustrated by the inset, where the black line shows the synaptic activation function and the horizontal red lines give the position and response range of V _s in response to the square wave. When Δ V _s becomes more negative, V _s shifts to shallower parts of the activation curve. As a result, the gain decreases, and thereby, the effectiveness of the feedback decreases as well. Therefore, the voltage swing of V _s gradually becomes larger (i.e., closer to that of V _p, the signal at the cone pedicle). The negative part of V _s leads to a strongly compressed negative part of V _h. This reduces or even abolishes the slow sagging, which is, in the model, attributed to local processes (in particular a nonlinear membrane) in the cone inner segment, axon, and pedicle ( Figure 1A; van Hateren, 2005). The positive part of V_s is still in a steeper part of the activation curve, and therefore, the positive part of V_h is left mostly intact, including the slow sagging occurring at off-steps. Only at strongly negative ΔV_s are both negative and positive parts of the response nearly vanished in V_h. 

Discussion

In this article, the previous temporal model of the cone–horizontal cell circuit in the macaque (van Hateren, 2005) has been extended to include spatial processing. The extended model is well able to describe all available measurements on the spatiotemporal processing in H1 cells. It provides quantitative insight into the origin of the main aspects of signal processing by H1 cells, particularly the gain controls, as well as the spatial and temporal filtering characteristics. Moreover, the model is well suited as a preprocessing module for analyzing downstream neurons. 

Three new processes had to be added to accomplish a good performance of the model. First, the local feedback gain from horizontal cell to each triad in the cone pedicle was made adaptive, depending on the local synaptic activation in each triad. Second, the spatial integration of the network of H1 dendrites was modeled as a superposition of short- and long-scale spatial filters, in accordance with Packer and Dacey (2002). Third, an adaptive low-pass filter, positioned at the cone pedicle/synapse, was added. 

The spatiotemporal model is more complex than the previous temporal one but is still very much simplified compared with the real physiological system. Because the elements added are slow and essentially fixed for stimuli of fixed field size, the present model is fully consistent with the previous temporal model, and I verified that it can indeed be fitted equally well to all measurements used in van Hateren (2005). Nevertheless, the added complexity of the adaptive elements considerably complicates the interpretation of seemingly simple experiments addressing spatial integration by horizontal cells (cf. turtle: Benda, Bock, Rujan, & Ammermüller, 2001; cat: Lankheet, Frens, & van de Grind, 1990). 

Several simplifications were deliberately made during the development of the model. One simplification is that the cones are not electrically coupled in the model, whereas it is known that macaque cones are coupled in reality (Hornstein, Verweij, & Schnapf, 2004). I found that adding a realistic cone coupling to the model did not significantly change the model responses to the stimuli considered here, and I therefore omitted it from the model. However, it is likely that it would have to be added when, for example, the influence of noise on the processing is investigated (DeVries, Qi, Smith, Makous, & Sterling, 2002). A second simplification is that the spatial integration is modeled as separable from the temporal processing, whereas these are nonseparable in real two-dimensional cables. Furthermore, the simple superposition of short- and long-scale spatial filters seems a poor physical model of the electrical properties of coupled H1 cells, although simulations with a compartmental model by Packer and Dacey (2005) suggest that it may lead to phenomenologically realistic results. 

Although the primary aim of the present modeling is to obtain a model of the cone–H1 circuit that is working and is as simple as possible, a secondary aim is to have a physiological interpretation of all components. Whereas the previous temporal model could rely to a considerable extent on the detailed picture of phototransduction that has been drawn over the past decades, the processes in the cone triad have not yet been established with comparable confidence. Several mechanisms added to the model therefore have a less clear physiological interpretation, although they all appear to be well within the range of possibilities consistent with the present knowledge of the system. 

The most important addition to the model is an adaptive feedback gain. This mechanism has not yet been directly observed experimentally and should, therefore, be viewed as a hypothesis, that is, as a prediction of the model. Its effects should be clearly observable with stimuli that change the working point of the synapse by abrupt changes in the spatial or chromatic composition of the stimulus (see also below). Although the proposed mechanism appears to be the simplest explanation of the experimental data, it cannot be ruled out that it is a stand-in for a more complicated synaptic mechanism. Obviously, only targeted experiments can clarify this issue. 

The adaptive feedback gain was scaled between 0 and 1 by Equation 2 to be consistent with an ephaptic (extracellularly electrical) mechanism of feedback (Kamermans et al., 2001; Kamermans & Fahrenfort, 2004). However, it should be noted that the maximum gain can be made smaller or larger than 1, which might occur with a nonephaptic feedback mechanism, for example, through the H⁺ concentration (Vessey et al., 2005). Any change in the maximum gain of g_h can be formally compensated in the model by a reciprocal change in the synaptic gain g_t. This only scales V_p, V_d, and V_s (which have not been measured and are, therefore, undetermined to some extent) while keeping the measured V_h the same. 

Given a fixed maximum gain of g _h, the actual g _h is varied by the model depending on how much transmitter is released in the triad as a result of the particular stimulus presented. Assuming ephaptic feedback, such a change in g _h might be produced by a modulation of the hemichannels (e.g., by protons; Sáez, Retamal, Basilio, Bukauskas, & Bennett, 2005), a shift in other conductances of the horizontal cell dendrites, or a shift in conductances of the surrounding neurons and structures involved in the triad (e.g., through a slow GABAergic system; Kamermans, Fahrenfort, & Sjoerdsma, 2002; reviews: Kamermans & Spekreijse, 1999; Schwartz, 2002; Wu, 1992; but see Verweij et al., 2003). With a more conventional chemical feedback mechanism, there is a range of possibilities for gain control. However, the model assumes that the feedback is primarily subtractive, and I find it difficult to see how a strongly nonsubtractive feedback system (i.e., multiplicative or divisive) could be equally consistent with the measurements considered here. This implies that a chemical feedback mechanism needs to be approximately subtractive over the voltage range normally occurring in the feedback. Although this is certainly feasible, it constrains the ways by which such a mechanism can be realized. Note that ephaptic feedback is inherently subtractive and, therefore, a particularly interesting possibility for the physiological substrate of the model proposed here. 

A second mechanism added to the model, an adaptive low-pass filter, has its physiological substrate presumably in the cone pedicle or in an early part of the synaptic mechanism at the triad. In principle, a change in the membrane time constant of the cone pedicle might be a possibility but that is not likely here. Such a change would normally also change the presynaptic gain, and I found that linking the change in τ _p with such a presynaptic change in gain could not provide acceptable fits to the measurements. Therefore, either a more complicated process in the pedicle (possibly driven by the GABAergic system) or an electrically silent synaptic process restricted to each triad may be involved. A constraint is that it should be positioned before the subtraction ( V _p − V _d) because I found that assuming an adaptive time constant as part of the main feedback loop did not produce adequate fits. 

Afterimages and chromatic adaptation

The most conspicuous new mechanism proposed here, the adaptive feedback gain, is likely to have important functional consequences for visual perception. Before discussing that, it should be noted that it can be interpreted more mundanely as a (slow) mechanism that tries to keep the synapse close to its normal working range: Whenever the local stimulus is too bright, and therefore, V _s is too negative and the transmitter output of the triad is low, the increase in g _h will make V _d more negative and V _s more positive, mostly restoring normal transmitter output. Conversely, when the local stimulus is too dark, the high transmitter output will gradually decrease g _h and thereby make V _d more positive and V _s more negative, thus reducing transmitter output. Note that the (slow) change in g _h only affects the local triad and, therefore, assures that feedback remains effective despite large variations in loop gain that might otherwise result from the nonlinearity of the synaptic activation, I _t( V _s). Although such variations do occur on a fairly short time scale, steady variations are partly compensated by the change in g _h. 

Because the change in g _h is slow, it forms a sort of spatial visual memory of the stimulus, and therefore, it will contribute to the occurrence of negative afterimages. When the eye fixates a bright spot for a while, the local gain, g _h, increases. When, subsequently, a homogeneous gray field is viewed, the large g _h at the spot position will drive V _d more negative and, therefore, V _s more positive. Thus, I _t is larger at that position than in the surround and, therefore, perceived as dark, that is, as the dark afterimage of a bright spot. Because the adaptation of g _h is assumed to take place locally at the cone triad synapse, the afterimage retains the full spatial resolution of the cone sampling grid. Note that the dynamics of afterimage formation and fading are determined by τ _itd in the model. In reality, this process might have more complicated dynamics than a first-order low-pass filter, but it may prove difficult to model this in more detail as long as the underlying physiological processes of this mechanism are not identified. 

The adaptive feedback gain is also likely to contribute to chromatic adaptation. The H1 cell receives input from all L (red) and M (green) cones within its receptive field (Dacey, 1999). If a scene is reddish because of a reddish illuminant, the L cones will, on average (over time), be more strongly stimulated (i.e., more hyperpolarized) than the M cones. Because the H1 cell averages over L and M cones, I_t at the L cones will then be smaller than normal and I_t at the M cones will then be larger. This will gradually make g_h at the L cones larger (i.e., it will make V_d more negative and, thus, V_s more positive; hence, I_t will be larger and, thus, closer to normal) and g_h at the M cones smaller (hence, I_t will be smaller and, thus, closer to normal). The result is an I_t for both the L and M cones that is rather similar to the I_t that would arise from illuminating the scene with a neutral gray. In other words, the adaptation of g_h functions as chromatic adaptation, and the result is that perception (i.e., I_t) is partly compensated for shifts in the chromaticity of the illuminant. Chromatic adaptation with a time constant in the order of 10–20 s has, indeed, been observed in both psychophysical measurements (Rinner & Gegenfurtner, 2000) and macaque retinal ganglion cells (Yeh, Lee, & Kremers, 1996). 

Acknowledgments

I wish to thank Herman Snippe for numerous discussions and Barry Lee, Herman Snippe, Maarten Kamermans, and Duco Endeman for helpful comments on the manuscript. 

Commercial relationships: none. 

Corresponding author: Hans van Hateren. 

Email: [e-mail]. 

Address: Department of Neurobiophysics, University of Groningen, Nijenborgh 4, NL-9747 AG Groningen, The Netherlands.