The light intensity of visual stimuli varies enormously in natural environments. Visual systems generally cope with these variations by regulating sensitivity to light at an early stage of visual processing. In the visual system of the macaque, Smith, Pokorny, Lee, and Dacey (2001) showed that such a regulation could be measured at the level of the horizontal cells in the outer retina. These cells are directly postsynaptic to the cones. They are thought to provide feedback to the cone pedicles, and thus regulate the temporal, spatial, and spectral properties of the output of the cones to the bipolar cells.

A fully adequate model of sensitivity regulation in the outer retina is lacking, whereas it is clearly important for modeling visual processing upstream in the retina and beyond. It should necessarily form the first stage of any model of the visual system that could function well in natural, outdoor-lighting conditions. The present study, therefore, aims to develop such a model, using the extensive set of horizontal cell measurements recently published (Lee, Dacey, Smith, & Pokorny 2003; Smith et al., 2001). It will be shown that the sensitivity regulation observed in horizontal cells with wide field, spectrally white stimuli, is fully consistent with the regulation expected from the known processes of the phototransduction machinery of the cones.

Both phototransduction and the neural circuit of cone and horizontal cell have been extensively investigated, and much is known about the physiological processes that determine their dynamics (for reviews of phototransduction, see Arshavsky, Lamb, & Pugh 2002; Burns & Baylor, 2001; Fain, Matthews, Cornwall, & Koutalos 2001; Pugh & Lamb, 2000; for a review of horizontal cells, see Kamermans & Spekreijse, 1999). Below I will represent a series of these processes, first as differential equations describing their dynamics, and subsequently as a system model consisting entirely of first-order low-pass filters and static nonlinearities. Because many readers may not be fully familiar with all of these fields, I will present the key topics in a relatively self-contained manner. The variables in the model have different physical dimensions, ranging from trolands through ion currents and concentrations to electrical currents and voltages. To keep the equations as concise as possible, I will not explicitly include dimensional conversion constants, and I will allow gains that are not dimensionless. Wherever the dimensions in an equation appear not to match, the reader should assume an implicit conversion constant of unit dimension.

In the sections below I will show that the (coupled) differential equations describing the system can be translated into a system model with only two types of components: first-order low-pass filters and static nonlinearities. The fact that the behavior of these components is straightforward and well understood greatly simplifies understanding and analysis of the system dynamics. I will therefore first introduce these basic components.

The standard equation describing a low-pass filter with time constant τ_y, transforming an input signal x(t) into an output signal y(t), is

An alternative form of a low-pass filter that often occurs when dealing with chemical reactions and membrane processes is given by a variant of Equation 1:

The low-frequency gain now equals τ_y. An example of such a process, where gain and time constant are covarying, is the voltage response to a current injected into a cell: When the membrane resistance is lowered, both the membrane time constant and the gain (voltage response to a given current) decrease. Below, such processes are described as the cascade of a separate gain τ_y and a standard low-pass filter as in Figure 1A.

The second component used in the models below is the static (memory-less) nonlinearity depicted in Figure 1B: It gives, for each point in time, the output y(t) as a function f(x(t)) of the input x(t). The function f(·) itself is not a function of time, only of the present value of its input x. A special case is the function f(x) = gx, with g a constant describing a (linear) gain.

For the equations describing the phototransduction cascade I will keep the notation as short as possible by mostly using single letter symbols for the variables and constants involved. Furthermore, I will keep the number of variables, such as scaling constants, to a bare minimum, and include only those that affect the calculations when fitting to the measurements on horizontal cells. Cascaded scaling factors or gains will therefore be merged into a single factor or gain wherever possible. Many of the equations are based on modeling work by Calvert, Govardovskii, Arshavsky, and Makino (2002), Detwiler, Ramanathan, Sengupta, and Shraiman (2000), Hamer (2000), Koutalos and Yau (1996), Pugh and Lamb (2000), Schnapf, Nunn, Meister, and Baylor (1990), and in particular Nikonov, Engheta, and Pugh (1998), and Nikonov, Lamb, and Pugh (2000). These studies also give extensive references to the original experimental studies that have led to the understanding of the dynamics of the various parts of the system.

Figure 2A shows on the left the reaction scheme leading from the absorption of light (I) by the visual pigment (R) via a G-protein (G) to the production of activated phosphodiesterase (E^*). In the first step, the visual pigment absorbs light, and is converted into an active form (R^*). The active form is removed with a rate constant k_R. It is assumed here that the light intensity is sufficiently low such that the concentration of R is not significantly affected (new R is produced fast enough, and there is no significant bleaching of the pigment). The reaction is then described by

In the second step, R^* activates a G-protein into its active form, G^*, which subsequently forms an active complex E^* with phosphodiesterase (PDE). Although the latter reaction can be modeled as a separate low-pass filter, its time constant is assumed to be so short that it can be neglected (Lamb & Pugh, 1992). It is again assumed that light intensities are such that neither G nor PDE is significantly reduced. Therefore, the production of E^* is described by

Figure 2B shows on the left the reaction scheme where E^* reduces the concentration of cGMP (symbolized by X below), which is itself produced at a rate α. The activity of PDE, β, is described by a constant dark activity (here c_β, often called β_dark in the literature) plus a term depending on the amount of activated PDE:

Figure 3 shows on the left the reaction scheme of the calcium feedback loop. An increase in cGMP concentration (X) will open more cyclic nucleotide-gated (CNG) ion channels in the plasma membrane of the outer segment, resulting in an inward current consisting partly of Ca²⁺. The resulting increase in calcium concentration is counteracted by a Na⁺/Ca²⁺–K⁺ exchanger, and possibly by calcium buffering. The increased calcium concentration reduces the rate with which guanylate cyclase (GC) synthesizes cGMP, and thus counteracts the initial rise of cGMP. In a secondary loop, more important in cones than in rods (Rebrik & Korenbrot, 2004), Ca²⁺ decreases the sensitivity of the CNG channels to cGMP.

The Ca²⁺ concentration, called C here, is described by

The reaction scheme contains two nonlinearities due to cooperativity. The opening of the CNG channels depends on the cGMP concentration as

Equations 12-15 are represented by the system model on the right side of Figure 3. To emphasize that it is a negative feedback loop and to enable an easy comparison with existing models, it is drawn as a divisive gain control: The input Q (=1/β) is divided by 1/α, which is itself an expansive function of the calcium concentration, and therefore provides strong negative feedback. The broken line in the Figure represents the reduction of the channel sensitivity caused by calcium: either through the calcium concentration C or, possibly more accurately, through the calcium current if the interactions are restricted to local domains. To keep the model as simple as possible, I am not modeling this interaction in detail. Instead, I assume it accounts for the fact that the standard value n_X =3 cannot produce acceptable fits of the model to the measurements. It produces major deficiencies in the sensitivities as a function of contrast and mean intensity. Leaving n_X as a free parameter in the fits leads to much lower values, typically n_X =1. This low value of n_X may be the result of the direct calcium feedback onto the CNG channels: Channel openings caused by an increase in X are counteracted by the subsequent increase in calcium influx. It is well known in engineering that negative feedback loops can apparently linearize nonlinearities in the forward path. For example, a fast divisive feedback loop with a squaring operation in the feedback path will have output = input/output², or output = input^1/3, thus converting an actual n_X =3 into an apparent n_X =1. I will assume the value n_X =1 below.

Similarly, I found that the fits, although not unreasonable for values n_C = 2 or 3, were significantly better for n_C ≈ 4. An apparent value of n_C = 4 was recently reported and discussed for rods (Burns, Mendez, Chen, & Baylor, 2002). I will assume n_C =4 below.

The inner segment contains a range of ion channels, partly voltage sensitive, shaping the response (Yagi & MacLeish, 1994). To model this without introducing many (poorly known) parameters, I will follow here the approach of Detwiler, Hodgkin, and McNaughton (1980). It is assumed that the overall steady-state conductance of the membrane, g_is, is given by a (nonlinear) function of the receptor potential V_is:

The change of V_is in response to a change in I_os is in reality not as instantaneous as suggested by Equation 17, because the membrane capacitance has to be charged. This can be represented by an additional low-pass filter with time constant τ_m, assumed to be approximately constant. The gain of this filter is merged into the scaling constant a_is. The above equations lead to the system model shown in Figure 4. Although the signal transfer from the cone inner segment to the cone pedicle may produce some additional filtering (Hsu, Tsukamoto, Smith, & Sterling, 1998), it is assumed here that the signal arriving at the cone pedicle is that of the inner segment. An alternative interpretation is that any additional filtering can be thought to be incorporated into the parameters γ, a_is, τ_is, and τ_m used here for the model of the inner segment.

For the feedback from the horizontal cell to the cone pedicle, I assume the subtractive scheme shown in Figure 5A. This is consistent with recent proposals for the feedback mechanism (Kamermans, Fahrenfort, Schultz, Janssen-Bienhold, Sjoerdsma, et al., 2001), and also with the spatially local character of the gain control (see Discussion). The driving force V_s for transmitter release I_t is determined by the difference of the voltage of the inner segment, V_is, and the voltage V_h of the horizontal cell multiplied by a gain g_h. Thus

The release of transmitter and its postsynaptic transduction are known to be nonlinear processes. The cone-horizontal cell feedback loop, possibly together with more local ionic feedback circuits, will tend to (apparently) linearize this process (Kraaij, Spekreijse, & Kamermans, 2000). The model as presented up to this point, with a linear gain g_s, produces acceptable fits to all horizontal cell measurements considered below. There are two nonlinearities that improved the fits sufficiently to justify their inclusion into the model. For high-contrast steps at high illuminance, there is an indication of a transiently reduced gain g_s (see Comparison with measurements, Figure 8). To model this, I included a nonlinearity similar to the one used by Kamermans, Kraaij, and Spekreijse (2001): The effective transmitter release as a function of V_s is, for the values of V_s relevant for the fits, described by a Boltzman function:

A second nonlinearity becomes clear when inspecting the high-frequency oscillations seen in response to pulses (Figure 7). The frequency of these oscillations decreases with decreasing background intensity. This indicates that the time constants or gain of the feedback loop are intensity dependent. This also influences the high-frequency part of the curves giving the sensitivity as a function of illuminance and frequency (Figure 14). I found that the shift in oscillation frequency could be modeled by assuming that a single factor, a_I, decreases the gain g_t and increases the time constants τ₂ (or, equivalently, τ₁) and τ_h with decreasing intensity. The effective gain and time constants are then g_t/a_I, a_Iτ₂, and a_Iτ_h. Decreasing intensity covaries with increasing V_is; Good fits were obtained by making a_I depend on V′_is, a low-pass filtered version of V_is, as

The model is shown in its complete form in Figure 6A. Figure 6B shows how a simple stimulus, a 100-ms step of contrast 2 at a background illuminance of 100 td, is transformed by the subsequent processing steps in the model, using the parameters fitted to Figure 7 below. The encircled numbers in the panels correspond to the numbers in Figure 6A. The intensity step is first low-pass filtered by τ_R and τ_E, yielding a signal proportional to the concentration of activated PDE (panel 2). The rate of cGMP hydrolysis, β, is elevated by c_β relative to the E^* curve (Equation 10), hardly visible in the graph at these illuminances (panel 3). The inverse of β (panel 4) is the signal that acts as the input to the calcium feedback loop. It is regulated by 1/α (panel 8), a low-pass filtered version of the loop output. This signal is slightly delayed relative to 1/β because of the low-pass filtering, and therefore boosts high temporal frequencies (panel 5). Furthermore, the feedback effectively reduces the dynamic range needed by the signal (cf. panels 6 and 4). The feedback loop of the inner segment (panel 10) reduces low frequencies (panel 9). In panel 9 the membrane voltage of the cone is shown relative to its membrane potential in the dark (i.e., V_is-V_is,dark). Similarly, panel 14 shows the membrane voltage of the horizontal cell relative to its dark value (i.e., V_h-V_h,dark). It should be noted that the scaling of V_is relative to V_h is not fully determined by the measurements on V_h considered here, because it depends on fixing g_h = 1 in Equation 20. In the cone-horizontal cell feedback loop (panels 11-14), the membrane potential of the horizontal cell is subtracted from that of the cone. The resulting V_s = V_is-V_h is small, only in the order of a mV. It is transient, because the signal in the horizontal cell is delayed, due to the low-pass filtering, relative to that in the cone. The resulting oscillations (panels 11, 12) arrive, however, strongly attenuated in the horizontal cell (panel 14). Note that the signal presumably going to the bipolar cells (bc in Figure 6A) is considerably more transient (panel 13) than that of the horizontal cell. How much of this shows up in the bipolar cells depends on the amount of low-pass filtering and further processing occurring in the bipolar cells.

I implemented the model as a series of auto-regressive moving average (ARMA) filters (e.g., Fante, 1988; also see below). A Fortran90 code producing the example of Figure 6B can be obtained by clicking here (Supplementary material). I verified the correctness of the model implementation by comparing its responses to small signals with the transfer function obtained from a small-signal analysis of the model (see Supplementary material), and by comparing its responses to large signals with a fourth-order Runge-Kutta solution of the system of differential equations involved. The implementation using ARMA filters is advantageous because it is both highly accurate and very fast, and therefore allows extensive fitting of the model to the horizontal cell data. The fitting was performed using a simplex algorithm (Press, Teukolsky, Vetterling, & Flannery, 1992) for minimizing the RMS-deviation between model responses and measurements. When responses with very different amplitudes r_max were used for simultaneous fitting, the RMS-deviations were scaled by to prevent the largest responses to completely dominate the fit. This was a pragmatic choice that still assigns more weight to large responses because they are less affected by noise and provide more information on the nonlinearities in the model. The measured data from horizontal cells in Figures 7-14 were obtained from the graphs in the on-line versions of Smith et al. (2001) and Lee et al. (2003), using specialized software (g3data).

For a first-order low-pass filter with time constant τ, the output y(n) to an input x(n), with a time step Δt, is given by the following ARMA filter (Brown, 2000):

Below I present fits of the model to a range of measurements on macaque H1 horizontal cells, as published in Lee et al. (2003) and Smith et al. (2001). The behavior of the model is complex, due to the nonlinearities and feedback loops, and the model contains many parameters. There is always the possibility, therefore, that good fits are obtained by overfitting (i.e., with very different, nonphysiological parameter values for different stimuli or different cells). To show that this is not a significant problem here, I determined a set of typical parameter values (listed in Table 1), which together determine a "generic model" producing generic responses. For most measurements presented below (black symbols), I will show two model calculations: the generic response (dashed blue line) and the fitted response (continuous red line). I often fitted with only a subset of parameters left free, fixing the others because which parameters are relevant depends on the type of stimulus. For example, the parameters determining the behavior of the model for large changes of background intensity should not be left free to vary for a medium contrast stimulus at one particular background intensity. Details on the fitted parameters and their values are given in the Figure captions and in Table 1.

Because the properties of horizontal cells, such as their sensitivity, vary, the generic response will generally deviate from the measurements. But as will become clear below, the generic responses always capture the basic qualifying characteristics of the measured responses.

Apart from the parameters of the model of Figure 6A, there is one additional parameter that I used as a free parameter in the fits where applicable, namely an overall delay. This accounts for any delay not incorporated into the model (such as transmission and diffusion times and possibly instrumental delay). The delay was forced to be identical for the responses to all stimulus conditions used for a particular cell. It was always small, typically around 3 ms.

The light intensities in this study are expressed in units of trolands (td). This is a measure of retinal illuminance, defined as the area, in mm², of the pupil of the eye times the scene luminance, expressed in candela/m². At a luminance of 100 cd/m², the pupil has a diameter of approximately 3 mm, thus then 100 cd/m² corresponds to 640 td. Such a luminance corresponds roughly to the mean luminance outdoors on a dull cloudy day, and is an order of magnitude higher than the typical mean luminance indoors.

Figure 7 shows responses to 10-ms pulses and 100-ms steps, with various contrasts superimposed on 100-, 10-, and 1-td backgrounds. Note the nonlinearity of response size as a function of contrast and of background intensity. The model fits (continuous red lines) were made to all stimulus conditions simultaneously, and are generally good, considering the wide range of stimulus conditions. The generic model responses (dashed blue lines) capture the qualitative characteristics of the measured responses.

It is possible to identify the parts of the model responsible for various aspects of the response shapes. The response size as a function of contrast and background intensity is mainly determined by the model components up to the generation of the photocurrent, I_os. The response sagging during the step responses at 100 td, and the response rebounds after pulses and steps, are mainly due to the properties of the model components representing the inner segment. The high-frequency oscillations are due to the cone-horizontal cell feedback loop.

At 100 td these oscillations are less prominent at high-contrast steps than at low-contrast steps, which points to the existence of a nonlinearity in the cone-horizontal cell feedback loop. This was modeled in Figure 5B. Figure 8 shows the effect of the absence of this nonlinearity (as in the model of Figure 5A). The effect of the nonlinearity can be understood as follows. During a high-contrast step, V_s becomes strongly negative, and thus brings the transmitter release, Equation 21, into a part of the curve with decreased slope. This decreases the small-signal gain, and this decreased gain in the feedback loop leads to a less steep rise of the response and less prominent oscillations. The fact that the oscillations in the model response are not as strongly reduced as in the measured response suggests that there are additional nonlinearities present, possibly due to voltage-sensitive channels in the horizontal cell membrane active at large hyperpolarizations.

I found that the best fits to the stimuli here and below were obtained with a low value of the time constant of the calcium feedback loop, with τ_C typically 3 ms. Forcing τ_C to be much larger worsened the fits considerably, and led to strongly biphasic cone photocurrents and membrane voltages in response to pulses and steps. This biphasic behavior is not consistent with the horizontal cell measurements considered here, which show only a mild rebound (Figure 7), with dynamics differing from those predicted from a slow calcium feedback. In rods, the photocurrent is generally not biphasic, unless the calcium dynamics are manipulated (Torre, Matthews, & Lamb, 1986). Although there are reports in the literature of strongly biphasic photocurrents and membrane voltages in primate cones (Schnapf, Nunn, Meister, & Baylor, 1990; Schneeweis & Schnapf, 1999), these findings may well be a consequence of a disturbed calcium dynamics due to the experimental techniques used. Such a disturbance of the cones is unlikely to have occurred in the horizontal cell measurements considered here, because cones were not directly manipulated, and the preparation left the retina mostly intact (Dacey, 1999; Smith et al., 2001).

Figure 9 shows responses to sinusoids of various frequencies at contrasts 0.25, 0.5, and 1, all at a 1000-td background level. In particular at high contrasts there are several distortions visible in the experimental data, which are also produced by the model. The compressive/expansive distortion at 0.61 Hz is mainly produced by the static nonlinearity, 1/β. The low-intensity part of the sinusoidal stimulus produces a small β, which is subsequently blown up by 1/β to the high peak in the response. The peak height is limited by the minimum value of β, c_β in Equation 10. The high-intensity parts of the sinusoid produce a large β, which is subsequently compressed by 1/β, which acts then as a compressive nonlinearity. The detailed shape of the distortion at 0.61 Hz is also determined by the calcium feedback loop, which in effect relinearizes the response to some extent: the high peak (low intensity, large 1/β) produces high levels of calcium, reducing the gain of the forward path (divisive gain control in Figure 3). This brings the response considerably closer to the steady-state level (dashed line). On the other hand, the low response (due to high intensity, small 1/β) produces low levels of calcium, increasing the gain of the forward path, also resulting in a response closer to the steady state than would have resulted without the calcium feedback. Because the former effect is stronger than the latter, the distortion is reduced.

For higher frequencies, the distortion due to 1/β gradually disappears, because the first two low-pass filters, τ_R and τ_E, reduce the depth of modulation of β and thus 1/β. At 4.88 and 9.76 Hz, another distortion becomes clearly visible: The falling flank of the response becomes steeper than the rising flank. This distortion is mainly due to the low-pass filters in the calcium feedback loop: The maximum reduction due to the control signal 1/α (Figure 6A) is only reached right after the peak in 1/β, which results in a steep falling flank right after that peak.

Finally, a third distortion can be seen at 30.3 Hz, where the rising flank is steeper than the falling flank. This is due to the filtering by the cone-horizontal cell feedback circuit. This circuit has a resonance frequency around 30-40 Hz, with strong phase changes near the resonance peak. This strongly shifts the phase of the distortion products, harmonics at multiples of the fundamental frequency, already produced by the cone. The first harmonic thus gradually shifts phase relative to the fundamental when going through the frequency range surrounding the resonance frequency, which happens to result in a steepening of the rising flank at 30.3 Hz. The distortions described here appear not to be specific to macaque H1 cells, because very similar distortions were measured in cat horizontal cells (Lankheet, van Wezel, & van de Grind, 1991).

It may be noted that some of the remaining deviations between measurements and fits in Figure 9 are related to small vertical offsets between them. Indeed, the fits improve when allowing for small errors, in the order of a mV or less, in the estimates of the steady-state potentials (horizontal dashed lines). Such errors could arise, for example, by small drifts in the intracellularly recorded membrane potential during the experiment.

In Lee et al. (2003) an experiment was performed to test the speed of the sensitivity regulation. It consists of a high-frequency test sinusoid superimposed on a low-frequency vehicle wave of high modulation depth. The local response to the test modulation then gives a measure of the sensitivity at a particular phase, and therefore illuminance, in the vehicle wave. Figure 10 shows measurements and model responses. As can be seen, the sensitivity regulation is almost instantaneous, with only a small delay mainly due to the low-pass filters in the calcium control loop. The test response can be quantified by extracting the amplitude of the fundamental frequency for each response cycle (Lee et al., 2003). Figure 11 shows for two cells how this response varies as a function of vehicle contrast. Clearly, the responses are well captured by the model.