Abstract

A computational model of human cones for intensities ranging from 1 td up to full bleaching levels is presented. The model conforms well with measurements made in primate horizontal cells, follows Weber's law at high intensities, and performs range compression consistent with what is known of cones in other vertebrates. The model consists entirely of processes with a clear physiological interpretation: pigment bleaching, saturation of cGMP hydrolysis, calcium feedback on cGMP synthesis, and a nonlinear membrane. The model is implemented according to a very fast computational scheme useful for simulations, and sample programs in Matlab and Fortran are provided as supplementary material.

Introduction

Human cones operate over a large range of luminances, ranging from mesopic (10

^{−3}to 10 cd/m^{2}) to photopic (10 to >10^{6}cd/m^{2}), corresponding to natural light levels ranging from starlight to sunlight (Hood & Finkelstein, 1986, Table 5.2). They accomplish this by a combination of slow and fast gain controls, ensuring a good response to even small contrasts at any background luminance. In this article, we will present a computational model of human long (L)- and medium (M)-wavelength cones that covers a considerable part of this luminance range, namely, from mid-mesopic (1 td, approximately white paper illuminated by moonlight) to high photopic (>10^{6}td, brighter than white paper illuminated by sunlight). The troland (td) is a measure of retinal illuminance, defined as the luminance in candelas per square meter multiplied by the pupil area in square millimeters. Below, we will mainly use the term “(light) intensity” when referring to retinal illuminance.The cone model is based on an earlier version developed for describing macaque horizontal cells in the 1- to 1,000-td range (van Hateren, 2005) and is extended here primarily by adding the dynamics of pigment bleaching as recently described in detail by Lamb and Pugh (2004) and Mahroo and Lamb (2004). The model is fully based on molecular and cellular processes known to occur in cones, in particular in the phototransduction machinery in its outer segment. We will show that the predicted responses of the model conform well with what is expected of cones and that the model produces several of their key features: an integration time that depends on light intensity, dynamic range compression, and, over much of the intensity range, a conformance to Weber's law (i.e., equal contrasts yield equal responses independent of adapting light intensity). The model is strictly deterministic and does not include intrinsic sources of noise in the cone.

A major reason for developing the present model is that it can be used as a general-purpose first module in models addressing visual information processing at downstream parts of the retina and beyond. To our knowledge, it is the first comprehensive model that spans the entire intensity range relevant for natural, daylight vision. The model may also be used for finding upper bounds on the role of gain control by the cones when studying higher visual processing. It is therefore put in a form that allows very fast computations, as needed for computing large cone arrays, and sample programs in Matlab and Fortran that perform these for the interested reader are provided.

Model

The model is an extension of the cone model developed for describing responses of macaque H1 horizontal cells between 1 and 1,000 td (van Hateren, 2005) and is congruent with other recent studies that model cones or rods (Hamer, 2000; Hamer, Nicholas, Tranchina, Liebman, & Lamb, 2003; Nikonov, Lamb, & Pugh, 2000; Tranchina, Sneyd, & Cadenas, 1991). The first part of the model was also successfully used for describing photocurrents in human cones (van Hateren & Lamb, 2006) as derived from the electroretinogram (Friedburg, Allen, Mason, & Lamb, 2004). Because the H1 cell is driven by the L- and M-cones, the model is only valid for those cones. It is unknown if and to what extent it applies to S-cones. For extending the model into the high-photopic luminance range, two additions proved to be necessary. The first is bleaching of cone pigment (Lamb & Pugh, 2004), and the second is a saturation of cGMP hydrolysis.

As in van Hateren (2005), boxes containing a “

*τ*” in Figure 1 represent first-order low-pass filters of unit DC gain and time constant*τ*. The gains of the physiological substrates of such filters are combined as much as possible into single parameters because this simplifies the analysis of the system (see the 1) and also makes it more apparent which parameters are important when fitting the model to measurements (van Hateren, 2005). Boxes in Figure 1 not containing a*τ*represent static linear or nonlinear operations on the variables as indicated.Figure 1

Figure 1

Pigment bleaching

In the first stage of the model shown in Figure 1A, cone pigment is excited by light and subsequently bleached. Recently, new measurements have indicated that the conventional description of bleaching recovery by first-order kinetics (Rushton & Henry, 1968) should be replaced by rate-limited kinetics (Lamb & Pugh, 2004; Mahroo & Lamb, 2004). We will follow here the theory as presented in Mahroo and Lamb (2004). In their formulation, the excitable pigment

*R*_{N}, the excited pigment*R*_{N}*, and the bleached pigment*B*are all normalized, that is, expressed as fractions of the total amount of pigment in the cone. This implies*R*_{N}+*R*_{N}* +*B*= 1. The sensitivity of the cone is then proportional to the fraction of excitable pigment,*R*_{N}= 1 −*B*−*R*_{N}*. The normalized*R*_{N}* follows from*R*_{N}* =*c*_{N}*R**, where*R** is a variable in the model (Figure 1A) proportional to the amount of excited pigment and*c*_{N}is a normalization constant (see the 1 for units and default parameter values).The intensity where the dot denotes time differentiation. Note that

*I*produces excited pigment*R** according to the following differential equation:$ \tau R R \u02d9*=I(1\u2212B\u2212 c NR*)\u2212R*,$

(1)

*R** is expressed here in the same units as*I*with unit gain for the*τ*_{ R}low-pass filter, where the gain and dimensional conversion are combined with other gains into the parameter*k*_{ β}( Figure 1A; see also van Hateren, 2005). A similar strategy is followed in many of the equations below.The bleached pigment is described by rate-limited kinetics (Mahroo & Lamb, 2004): The first term on the right-hand side represents the production of

$ B \u02d9= c NR*/ \tau R\u2212 K B \tau B , 0 B B + K B.$

(2)

*B*through the decay of*R**, and the second term represents the reconversion of*B*into excitable pigment, where*K*_{B}is a semisaturation constant and*τ*_{B,0}is the time constant of reconversion. If*B*≪*K*_{B}, the latter term is simply −*B/τ*_{B,0}, consistent with earlier first-order descriptions of pigment kinetics (Rushton & Henry, 1968). However, if*B*≫*K*_{B}, the term approaches*−K*_{B}/*τ*_{B,0}, which shows that the reconversion of*B*is then rate limited. In Mahroo and Lamb (2004), it is argued that this may be caused by a diffusion barrier for transporting 11-*cis*-retinoid from the retinal pigment epithelium to the cones, where 11-*cis*-retinal subsequently recombines with opsin. Equation 1 is depicted in the scheme of Figure 1A as the upper control loop, where the gain of This equation shows that the slow gain control provided by

*R** production is proportional to*R*_{N}= 1 −*B*−*c*_{N}*R**. This factor depends on the bleaching history, described by the lower control loop. This loop represents the dynamics of*B*by rewriting Equation 2 as$ \tau B B \u02d9 = g B R * \u2212 B , w \u2062 i \u2062 t \u2062 h \tau B = \tau B , 0 B + K B K B a \u2062 n \u2062 d g B = c N \tau B / \tau R .$

(3)

*B*has a variable time constant: For small*B**τ*_{ B}=*τ*_{ B,0}= 25 s (Mahroo & Lamb, 2004), whereas for large*B**τ*_{B}= 150 s (using*K*_{B}= 0.2 and*B*close to 1). Representing the pigment dynamics according to Equation 3 lends itself to a very fast computational scheme using autoregressive filters (see Model implementation). This scheme simulates the cone up to two orders of magnitude faster than numerically solving the differential equations by standard integration routines (see also the 1 and the supplementary material).Using the subscript “0” to denote the steady-state values of variables, we can find the steady state of Equation 1 from or The steady state of Equation 2 produces a quadratic equation in and consequently (using Equation 5, and This implies that the excited pigment

$0= I 0(1\u2212 B 0\u2212 c N R 0*)\u2212 R 0*$

(4)

$ R 0*=(1\u2212 B 0) I 0 1 + c N I 0.$

(5)

*B*_{0}when the*R*_{0}* of Equation 5 is substituted. The general solution (see the 1, Equation A22) has two interesting limits. The first limit occurs for*c*_{N}*I*_{0}≪*K*_{ B}*τ*_{ R}/*τ*_{ B,0}or*I*_{0}≪ 6,600 td. Then, it can be shown that*R*_{0}≈ 1 and that the “pigment bleaching” part of Figure 1A reduces to a single low-pass filter*τ*_{ R}; that is, the control loops are ineffective. The second limit occurs for*c*_{N}*I*_{0}≫*K*_{ B}*τ*_{ R}/*τ*_{ B,0}or*I*_{0}≫ 6,600 td, where Equation A22 reduces to$ B0=1\u2212 \tau R \tau B , 0 K B 1 + K B c N I 0 + 1 c N I 0,$

(6)

*R*_{N,0}* =*c*_{N}*R*_{0}*, and*R*_{N,0}= 1 −*B*_{0}−*R*_{N,0}*)$ R N , 0*= \tau R \tau B , 0 K B 1 + K B$

(7)

$ R N , 0= \tau R \tau B , 0 K B 1 + K B 1 c N I 0.$

(8)

*R*_{N,0}* goes to a fixed fraction of about 2.3 × 10^{−5}(∼1,100 molecules) at high intensities, whereas the available excitable pigment*R*_{N,0}becomes inversely proportional to*I*_{0}. Because*R*_{N,0}represents the initial gain of the phototransduction cascade, this inverse proportionality basically guarantees Weber's law at high luminances. Weber's law breaks down, however, at extremely low frequencies or at extremely high intensities (*I*_{0}≫ 10^{6}td; see Equation A13 in the 1 and Figure 5).Saturation of cGMP hydrolysis

Gain control at high intensities (>10,000 td) is dominated by pigment bleaching, and at low to medium intensities (1–1,000 td), it is accomplished by a combination of cGMP hydrolysis and calcium gain control (van Hateren, 2005). The latter two processes also explain the changes in the cone integration time observed between 1 and 1,000 td. Combining these processes with pigment bleaching, however, does not provide a good description of the cone response over the entire intensity range from 1 td up to full bleaching levels. In particular, we found that the combined processes were not able to satisfactorily bridge the gap between 1,000 and 10,000 td, where they fail to produce responses consistent with Weber's law. Although pigment bleaching could, in principle, accomplish Weber's law down to 1,000 td, this would imply that strong bleaching starts clearly earlier than has been established experimentally. Conversely, the cGMP hydrolysis and calcium gain control could not be extended to 10,000 td without compromising the required changes in gain and integration time between 1 and 1,000 td.

We found that the simplest way to solve this problem is to assume a modest additional gain control between the stages of pigment bleaching and cGMP hydrolysis. In Figure 1A, this process is tentatively represented as “saturating cGMP hydrolysis.” Saturation might occur when a large amount of where

*E** is produced, and cGMP cannot access*E** fast enough to keep all*E** hydrolyzing at full speed. This would then effectively limit the hydrolysis rate that can be obtained. This is represented in the model in a simplified form by a compressive nonlinearity, giving the effective hydrolysis rate*β*_{e}as a function of the hydrolysis rate*β*without saturation:$ \beta e= \beta 1 + \beta / \beta e , max,$

(9)

*β*_{e,max}is a constant. Note that*β*_{e}=*β*when*β*≪*β*_{e,max}, and then this part of the model is identical to the 1- to 1,000-td version described in van Hateren (2005). For*β*≫*β*_{e,max},*β*_{e}goes to a maximum value of*β*_{e,max}.The above mechanism of saturating cGMP hydrolysis is quite hypothetical and certainly not the only possibility. Figures 1B and 1C show alternatives. In Figure 1B, it is the production of excited PDE (

*E**) that is a limiting factor. This could occur because there is an excess of*R*over*E*(Pugh & Lamb, 2000), and therefore, a large amount of*R** may saturate the production of*E**. The details of such a reaction scheme are not known, but it is likely that it would produce a compressive nonlinearity (symbolized by NL in Figure 1B) that is similar in shape to that of Equation 9. Although this nonlinearity is at a slightly different position in the transduction chain than the one in Figure 1A, it would lead to nearly the same response characteristics. Similar response characteristics are also expected from the third scheme, as shown in Figure 1C. This scheme involves a regulation of cone pigment phosphorylation by calcium, affecting both the time constant*τ*_{R}and the associated gain*g*_{R}through a dynamic nonlinearity (DNL). Regulation of pigment phosphorylation is well known in rods and has been modeled in detail (Hamer, Nicholas, Tranchina, Lamb, & Jarvinen, 2005). It has been reported in cones as well (Kennedy, Dunn, & Hurley, 2004), but its role in primate cones has not yet been firmly established.Calcium feedback and inner segment

The remaining processes in the model cone of Figure 1A have been described before (van Hateren, 2005; van Hateren & Lamb, 2006; van Hateren & Snippe, 2006; see also Nikonov et al., 2000; Pugh & Lamb, 2000). Briefly, activated PDE,

*E**, hydrolyzes cGMP at a free-space rate*β*which consists of a dark rate,*c*_{β}, and a factor proportional to*E**. As argued above, the actual, saturating hydrolysis rate is*β*_{e}. Because the hydrolysis is governed by a nonlinear differential equation, the light level determines both the gain of the resulting filter (as 1/*β*_{e}) and its time constant (*τ*_{X}= 1/*β*_{e}). The 1/*β*_{e}factor can be shown to produce super-Weber behavior (van Hateren, 2005, supplementary material; see also van Hateren & Snippe, 2006), which is subsequently corrected by the calcium feedback loop (see below). The resulting current produced by the outer segment,*I*_{os}, is finally converted to a voltage over the inner segment membrane,*V*_{is}, with the final feedback loop in Figure 1A representing the nonlinearity of the membrane.Model implementation

All temporal filters in the model are first-order low-pass filters of the general form with where

*τ*$ y \u02d9$

= *x*−*y*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(n)= f 1y(n\u22121)+ f 2x(n\u22121)+ f 3x(n),$

(10)

$ f 1 = exp ( \u2212 1 / \tau \u2032 ) f 2 = \tau \u2032 \u2212 ( 1 + \tau \u2032 ) exp ( \u2212 1 / \tau \u2032 ) f 3 = 1 \u2212 \tau \u2032 + \tau \u2032 exp ( \u2212 1 / \tau \u2032 ) \tau \u2032 = \tau / \Delta \u2062 t ,$

(11)

*τ*is the time constant of the low-pass filter and Δ*t*is the time step (here taken as 0.1 ms). Fitting to the data in Figure 2 was performed using a simplex algorithm (Press, Teukolsky, Vetterling, & Flannery, 1992) for minimizing the RMS deviation between model responses and measurements.Figure 2

Figure 2

Details on the model equations can be found in the 1. The supplementary material provides Fortran and Matlab programs that implement the ARMA scheme, as well as an equivalent Matlab program using a standard solver for coupled ordinary differential equations, and those that implement a small-signal expansion of the model (see the 1). Results are independent of programming language and computational scheme.

Results

Weber's law in horizontal cells

The dark gray trace in Figure 2 is taken from Figure 5C of Lee, Dacey, Smith, and Pokorny (2003) and shows the average sensitivity of 15 macaque horizontal cells to a low-contrast test stimulus (a 19.5-Hz sinusoid) superimposed on slowly varying background intensities. The measurements encompass retinal illuminances between approximately 1 and 2,000 td. We extended this curve beyond 2,000 td by assuming that Weber's law applies to those intensities; that is, the extension (dashed blue line) has a slope of −1 and is connected to the end point of the measured data curve. The solid red line is a fit to the complete curve (measured and extension) with the H1 model of van Hateren (2005), enhanced by adding pigment bleaching and saturating cGMP hydrolysis. The fit thus provided an estimate of the parameter

*β*_{e,max}as needed in Equation 9 for describing saturation of cGMP hydrolysis. Although the fit was made with the H1 model, the model part that describes the horizontal cell and its feedback to the cones acts as a linear system for this stimulus (computed here as a low-contrast, 19.5-Hz test modulation superimposed on steady backgrounds of various intensities). Therefore, the shape of the curve, including the Weber behavior over much of its range, is fully determined by the properties of the cone model.For comparison, the model prediction for a test stimulus frequency of 0.61 Hz rather than 19.5 Hz is also shown. At low frequencies, the range of adapting intensities where the cone is displaying Weber behavior or near-Weber behavior (with a slope of about −0.7 rather than −1) is even larger than at 19.5 Hz.

Note that the curves at both frequencies show slight fluctuations in slope around the main trend. We believe that these fluctuations are unlikely to be observable in real cones because they are caused by the discrete (lumped) nature of the model in Figure 1A. In the real cone, many of the components are likely to become more diffuse by the fact that they represent distributed rather than lumped processes. For example, the calcium feedback may not be identical over the entire length of the cone because of differences in effective light intensity and differences in biophysical parameters. This would lead to slightly different fluctuations being contributed to the sensitivity curve by different parts of the cone and, therefore, a somewhat smoother curve.

Contributions to Weber's law

To obtain insight into which part of the cone model contributes at which background intensities to Weber's law, we show, in Figure 3, responses at different stages in the model of Figure 1A, all normalized to the response at 1 td. The stimulus was, again, a low-contrast, 19.5-Hz test stimulus at different adapting intensities. The black curve marked

*R** shows the response after the pigment bleaching stage in Figure 1A. At levels where pigment bleaching becomes significant, approximately >10^{4}td, the curve obtains a slope of −1 and, thus, conforms to Weber's law (see also Equations 8 and A13).Figure 3

Figure 3

Because

*E** and*β*are linearly related to*R** and because the curves are normalized, the curves for*E** and*β*(not shown) are identical to that of*R**. However, the red curve marked*β*_{e}shows that the saturation of cGMP hydrolysis changes the curve in the 10^{3}- to 10^{4}-td range. Although this change is rather modest, the compressive nonlinearity acting on*β*cannot be omitted without producing significant discrepancies with measured data, as argued above.Below the 10

^{3}- to 10^{4}-td range, the sensitivity is mainly determined by the concerted action of cGMP hydrolysis and calcium feedback. The dashed blue line shows the response after the 1/*β*_{e}nonlinearity in Figure 1A, where 1/*β*_{e}describes how the cGMP concentration depends on PDE activity. Note that 1/*β*_{e}is not directly observable in a normally functioning cone because it would be modified by the gain*α*of the calcium control loop before it leads to observables, such as the cGMP concentration (*X*in the model) and the membrane current*I*_{os}. It would be observable, though, in a cone where the calcium feedback has been made inoperative, for example, by clamping the Ca^{2+}concentration by buffering. As argued before (van Hateren, 2005; van Hateren & Snippe, 2006), the 1/*β*_{e}nonlinearity provides super-Weber behavior, leading to a slope of −2 rather than the required −1 (Figure 3). However, the gain provided by calcium feedback onto cGMP production by guanylylcyclase,*α*corrects this. The solid blue line marked*α*shows how this gain increases in the 10- to 1,000-td part of the curve. The resulting*I*_{os}(dashed red curve) is close to Weber's law over much of its range, which is also true for the membrane potential*V*_{is}.Adaptation of contrast responses

Adaptation of cone sensitivity has often been studied by measuring the response size to flashes of positive and negative contrast given on different adapting backgrounds (Burkhardt, 1994; Normann & Perlman, 1979). Figure 4 shows the performance of the cone model to this stimulus protocol. The inset shows a typical response, here the response to a 100-ms step to 200 td on a background of 100 td. The abscissa shows time in milliseconds, whereas the ordinate shows the response (

*V*_{is}) relative to the dark*V*_{is}, normalized by the response that would result when*I*_{os}= 0. Thus, the normalized response equals 0 in the dark and approaches −1 when nearly all transduction channels are closed at high light intensities. The main panel of Figure 4 plots the peaks of the responses to similar 100-ms flashes of a range of contrasts given on a range of backgrounds. To comply with convention in the literature, we inverted the sign of the response; that is, positive contrasts produce positive responses relative to the steady state, and negative contrasts produce negative-going responses. The red line in the figure shows the steady-state membrane potential (again, normalized and inverted). The black lines show the response peaks to negative and positive contrasts of up to 10^{2}or more, for background intensities of 1, 30, 100, 300, 1,000, 3,000, 10^{4}, 3 × 10^{4}, and 10^{5}td. The background intensity of each black curve coincides with the point where it crosses the red curve.Figure 4

Figure 4

Although corresponding curves have not yet been measured in primate cones or horizontal cells, and, therefore, the validity of the model to this type of stimulus cannot be directly evaluated, the curves in Figure 4 show a strong qualitative resemblance with curves measured in turtle cones (Burkhardt, 1994; Normann & Perlman, 1979) and with curves for primate cones obtained through an extracellular technique (Valeton & van Norren, 1983). At low background intensities, the response range is mainly used for positive contrasts, whereas this gradually changes toward high intensities. The steady state gradually rises, until it reaches a plateau around 10

^{4}td, caused by pigment bleaching and the saturation of cGMP hydrolysis. From that point on, the curves are merely shifting to the right when background intensity is increased. Note that even high contrasts do not reach a normalized response of 1. This is a consequence of the saturating cGMP hydrolysis, limiting the maximum value that*β*_{e}can reach. A small fraction of the transduction channels therefore remains open even at the highest contrasts. Also note that there is a clear asymmetry in the response curves at high intensities, favoring negative contrasts over positive contrasts.Notwithstanding the resemblance, several properties of turtle cones are not captured by the cone model. The steady-state response in turtle cones goes at high intensities to a value closer to 0.5 than in Figure 4. Furthermore, turtle cones appear to reduce their dynamic range somewhat at high luminances (Burkhardt, 1994), which is also not present in Figure 4. However, these differences are relatively minor compared with the overall correspondence between the turtle measurements and the present model calculations. Moreover, it is not clear at present whether these differences represent a deficiency of the model or its parameter settings or a genuine difference between human and turtle cones.

Frequency sensitivity

Using the small-signal analysis of the model (see the 1), the frequency sensitivity curves of the cone can be readily computed. Figure 5 shows the sensitivity at a range of background intensities (numbers at the curves, in trolands). Several properties that are well known in cones are represented by these curves. The gain (millivolts per troland) decreases with increasing intensity, and for high intensities, this occurs inversely proportional to the intensity, at least for midrange frequencies. As a result, the sensitivity to contrast will become independent of intensity (Weber's law) because a fixed contrast produces intensity modulations proportional to the background intensity. Figure 5 shows that this is valid over a larger range of intensities at midrange frequencies (gray area) than at high frequencies. This is related to the fact that the cutoff frequency of the cone increases with increasing intensity. Note that the high-frequency parts of the curves nearly coincide at low intensities (i.e., below intensities where bleaching or cGMP hydrolysis limits occur). This so-called high-frequency linearity is a well-known property of temporal processing by the human visual system (Kelly, 1961).

Figure 5

Figure 5

Finally, the figure also shows a less well known consequence of bleaching kinetics. At very low frequencies and at very high luminances, Weber's law breaks down, as can be seen most clearly as the low-frequency falloff in the 10

^{5}- and 10^{6}-td curves (see Equation A13 in the 1). The gradual drop at very low frequencies is related to the fact that at very high intensities, the steady-state amount of excited pigment becomes constant (∼1,100 molecules, see above, following Equation 7). This necessarily implies that the sensitivity of cones at very high intensities is 0 for a frequency of 0 Hz. Note, however, that this breakdown of Weber's law is of little consequence for human visual perception because for frequencies of importance to perception (such as the 0.2- to 30-Hz range marked by the gray area in the figure), Weber's law holds quite well up to the highest intensities shown.The curves presented in Figure 5 bear a strong resemblance to psychophysical sensitivity measurements in the 2.5- to 40-Hz range recently reported by Stockman, Langendörfer, Smithson, and Sharpe (2006). The remaining differences are probably related to the fact that the psychophysical measurements engage more of the visual system than just the cones. The stronger low-frequency roll-off in these measurements is most likely related to similar filtering known to occur in the bipolar cells (Dacey et al., 2000) and ganglion cells (Lee, Pokorny, Smith, Martin, & Valberg, 1990) of the magnocellular pathway. The stronger high-frequency falloff in the psychophysical measurements may originate in a cortical low-pass filter as postulated by Lee et al. (1990).

Discussion

In this article, we have presented a model for human L- and M-cones that produces credible responses at intensities ranging from mid-mesopic levels (1 td) up to high-photopic levels (>10

^{6}td). To our knowledge, this is the first comprehensive model that spans the entire range of luminances that is of interest for natural, daylight vision.The model consists entirely of processes that have a clear physiological interpretation and that are mostly well established. By a cascade of the nonlinear processes of pigment bleaching, saturating cGMP hydrolysis, calcium feedback, and a nonlinear membrane, the cone accomplishes several major goals useful for vision. First and foremost, it compresses an intensity range of more than six orders of magnitude into the limited span of its membrane potential while retaining a quite high contrast sensitivity at each adapting background ( Figure 4). Second, it obtains contrast constancy (Weber's law) or near contrast constancy over much of its intensity range ( Figure 3), thus simplifying the downstream processing and interpretation of visual information. Finally, it regulates its integration time in such a way that it reduces the effects of shot noise caused by low numbers of photons and molecules at low light intensities ( Figure 5).

Although the model has been fairly well validated for the 1- to 1,000-td range in van Hateren (2005) and is supported by data obtained from human electroretinograms (van Hateren & Lamb, 2006), it should be noted that no such direct validation was possible for the extended intensity range discussed here. Pigment bleaching has been added purely based on the literature (Mahroo & Lamb, 2004), and it is not known to what extent cGMP hydrolysis in human cones is saturating. Furthermore, the fit in Figure 2 was made to partly extrapolated data. Finally, the model is strictly deterministic and does not include intrinsic sources of noise in the cone, which may be important when the limits imposed on visual computation are considered. Therefore, the present model should be considered as a first-order result, where it may prove necessary to adjust the default values of parameters or make alterations to the various modules, once more direct data on human cones in the 1- to 10

^{6}-td range become available.The major goal of the present model is to establish a tractable computational scheme that can be used as a preprocessing module for studying and modeling visual information processing in downstream parts of the human retina and beyond. Therefore, considerable efforts have been made to ensure that the structure of the model allows fast computations and is amenable to a straightforward small-signal analysis ( 1 and supplementary material). We believe that the model is well suited for this goal because it appears to capture all properties of cones that are likely to be perceptually important.

Supplementary Materials

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Supplementary File - Supplementary File

Appendix A

Model Equations

The computations for this article were performed using an ARMA scheme ( Equations 10 and 11) for each first-order low-pass filter, literally following the system diagram of Figure 1A. Matlab and Fortran programs that implement this scheme are provided as supplementary material. Alternatively, the system can be described as a series of coupled ordinary differential equations, which can be solved using a standard ODE solver. A Matlab program that implements this scheme is also provided as supplementary material. The two approaches yield identical results, but the ARMA scheme is considerably faster than the ODE scheme. The ARMA scheme is therefore particularly suited for simulating large cone arrays, when computing time becomes an issue. For the purpose of reference, the equations for the ODE scheme are listed below, with the following variables (and dimensions; “−” denotes dimensionless): where The following parameters and default values are used:

*t*time (in milliseconds);*I*retinal illuminance (in trolands);*R**, quantity of R* (in trolands);*B*quantity of bleached pigment as a fraction of the total amount of pigment (−);*E**, quantity of E* (in trolands);*X*concentration of cGMP (−);*C*free concentration Ca^{2+}(−);*I*_{os}, ion current through the outer segment channels (−);*V*_{is}, cone membrane potential (in millivolts);*g*_{i}, membrane conductance (1/mV);*β*free-space rate constant of hydrolysis of cGMP by PDE (1/ms);*β*_{e}, saturating rate constant of hydrolysis of cGMP by PDE (1/ms);*α*guanylate cyclase rate (implicit 1/ms).$ d \u2062 R * d \u2062 t=[I(1\u2212B\u2212 c NR*)\u2212R*]/ \tau R$

(A1)

$ d \u2062 B d \u2062 t= c NR*/ \tau R\u2212 K B \tau B , 0 B B + K B$

(A2)

$ d \u2062 E * d \u2062 t=[R*\u2212E*]/ \tau E$

(A3)

$ d \u2062 X d \u2062 t=\alpha \u2212 \beta eX$

(A4)

$ d \u2062 C d \u2062 t=( I o \u2062 s\u2212C)/ \tau C$

(A5)

$ d V i \u2062 s d \u2062 t=( I o \u2062 s/ g i\u2212 V i \u2062 s)/ \tau m$

(A6)

$ d g i d \u2062 t=( a i \u2062 s V i \u2062 s \gamma \u2212 g i)/ \tau i \u2062 s$

(A7)

$\beta = c \beta + k \beta E*$

(A8)

$ \beta e= \beta 1 + \beta / \beta e , max$

(A9)

$\alpha =1/[1+( a CC ) n C]$

(A10)

$ I o \u2062 s= X n X$

(A11)

*c*_{N}=*τ*_{ R}*K*_{cone}/*R*_{tot}= 4.1 × 10^{−9}(td)^{−1}, with*K*_{cone}= 60 isom s^{−1}(td)^{−1}and*R*_{tot}= 5 × 10^{7}molecules (Kenkre, Moran, Lamb, & Mahroo, 2005);*τ*_{R}= 3.4 ms (time constant of*R**);*τ*_{B,0}= 25 s (time constant of bleaching recovery);*K*_{B}= 0.2 (semisaturation constant of bleaching recovery);*τ*_{E}= 8.7 ms (time constant of*E**);*c*_{β}= 2.8 × 10^{−3}(ms)^{−1}(dark PDE activity);*k*_{β}= 1.4 × 10^{−4}(ms)^{−1}(td)^{−1}(*E** dependence of PDE activity);*β*_{e,max}= 4 (ms)^{−1}(semisaturation constant of*β*_{e});*n*_{X}= 1 (apparent Hill coefficient of CNG channel activation);*n*_{C}= 4 (Hill coefficient of GC activation);*τ*_{C}= 3 ms (time constant of Ca^{2+}extrusion);*a*_{C}= 0.23 (scaling constant of GC activation);*τ*_{m}= 4 ms (membrane time constant);*γ*= 0.7 (parameter of membrane nonlinearity);*τ*_{is}= 90 ms (time constant of membrane nonlinearity);*a*_{is}= 2.9 × 10^{−2}(mV)^{−1−γ}(parameter of membrane nonlinearity).Most of the equations above have been discussed in detail in van Hateren (2005). In particular, the equations for the inner segment, axon, and cone pedicle (Equations A6 and A7) should not be considered as a detailed biophysical model but rather as a simplified dynamic systems model. They summarize the nonlinear membrane properties of the cone by a single nonlinearity, following the analysis in the Appendix of Detwiler, Hodgkin, and McNaughton (1980). The time constant

*τ*_{m}is assumed to be approximately constant, representing not only the time constant of the inner segment membrane but also any additional low-pass filtering attributable to axon and pedicle (for further discussion, see van Hateren, 2005).The equations above produce the cone response without incorporating the small overall delay observed in cones (approximately 1–3 ms; van Hateren, 2005; van Hateren & Lamb, 2006). This delay can be straightforwardly incorporated, if desired. For the small-signal equations below, this would involve an additional transfer function exp(

*iωδ*), where*δ*is the delay.Small-signal analysis

For small signals, for example, small sinusoidal modulations used for obtaining a transfer function, it is much faster and more accurate to compute the transfer function directly from the analytical small-signal solution rather than using the full model equations. Moreover, the transfer function can be used, via the Fourier transform, to compute the response to arbitrary time-varying signals, provided that the modulations around the background intensity are small. For most of the processes in the model cone, a small-signal analysis has already been presented in the supplementary material of van Hateren (2005). Below, this is summarized and extended with the small-signal equations for pigment bleaching and saturating cGMP hydrolysis. Matlab and Fortran programs that implement these equations are provided as supplementary material.

We will first perform a small-signal expansion of Equation A1, where we assume The lower-case symbols represent small perturbations around a constant, steady-state value denoted by the subscript zero. Using similar methods as in the supplementary material of van Hateren (2005) and transforming to the frequency domain (with transformed variables denoted by a tilde), we obtain where the second line of the equation follows from

$ I = I 0 + i R * = R 0 * + r * B = B 0 + b .$

(A12)

$ r \u02dc * i \u02dc = ( 1 \u2212 B 0 \u2212 c N R 0 * ) 1 + c N I 0 + c N I 0 K B 2 \tau R ( B 0 + K B ) 2 \tau B , 0 + i \u2062 \omega \tau R + i \u2062 \omega \tau R \u2248 R N , 0 1 + i \u2062 \omega \tau R ,$

(A13)

*c*_{N}*I*_{0}≪ 1, with the exception of very low*ω*(see Figure 5).*c*_{N}*I*_{0}≪ 1 is true for normal intensities because 1/*c*_{N}> 10^{8}td. At low intensities, the fraction of excitable pigment,*R*_{N,0}, is close to 1, and then Equation A13 represents a simple low-pass filter. At high intensities,*R*_{N,0}is inversely proportional to the light intensity (Equation 8), and then the second line of Equation A13 represents Weber's law followed by a low-pass filter.Using Equation A3 yields and a small-signal analysis of Equations A8 and A9 gives Combining the results above with Equations 31, 32, 37, 45, and 50 of the supplementary material of van Hateren (2005) then gives the transfer function of the human cone as with In Equation A18, the with Here, The Once as needed in Equation A19, and as needed in Equation A18.

$ E * = E 0 * + e * \beta e = \beta e , 0 + \Delta \beta e ,$

(A14)

$ e \u02dc * r \u02dc *= 1 1 + i \u2062 \omega \tau E,$

(A15)

$ \Delta \u2062 \beta \u02dc e e \u02dc *= k \beta \beta e , \max 2 ( \beta e , \max + \beta 0 ) 2.$

(A16)

$ H c \u2062 o \u2062 n \u2062 e(\omega )= H R(\omega ) H E(\omega ) H \beta (\omega ) H \beta , o \u2062 s(\omega ) H i \u2062 s(\omega ),$

(A17)

*H*_{R}given by the first line of Equation A13,*H*_{E}by Equation A15,*H*_{β}by Equation A16, and*H*_{β,os}and*H*_{is}by$ H \beta , o \u2062 s= \u2212 n X \beta e , 0 + \beta e , 0 n X n C c 1 / ( 1 + i \u2062 \omega \tau C ) + i \u2062 \omega $

(A18)

$ H i \u2062 s= ( 1 + i \u2062 \omega \tau i \u2062 s ) V i \u2062 s , 0 ( 1 + i \u2062 \omega \tau i \u2062 s ) ( 1 + i \u2062 \omega \tau m ) + \gamma .$

(A19)

*β*of van Hateren (2005) is replaced by the equivalent*β*_{e}here. The steady-state values of variables required in these equations are obtained by numerically solving the following equation for*I*_{os,0}(or, equivalently, for*C*_{0}=*I*_{os,0}):$ I o \u2062 s , 0 1 / n X+ a C n C I o \u2062 s , 0 n C / n X=1/ \beta e , 0,$

(A20)

$ \beta e , 0 = \beta 0 1 + \beta 0 / \beta e , max \beta 0 = c \beta + k \beta E 0 * E 0 * = R 0 * R 0 * = ( 1 \u2212 B 0 ) I 0 1 + c N I 0 .$

(A21)

*I*_{0}is the steady-state intensity (in trolands).*B*_{0}follows from eliminating*R*_{0}* from the steady state of Equations A1 and A2, which yields$ B 0 = [ 1 \u2212 K B \u2212 \tau R K B ( c N I 0 + 1 ) \tau B , 0 c N I 0 + ( ( 1 \u2212 K B \u2212 \tau R K B ( c N I 0 + 1 ) \tau B , 0 c N I 0 ) 2 + 4 K B ) 1 2 ] / 2 .$

(A22)

*R*_{N,0}needed in Equation A13 then follows from$ R N , 0=1\u2212 B 0\u2212 c N R 0*.$

(A23)

*I*_{os,0}is obtained from Equation A20, one finds$ V i \u2062 s , 0=( I o \u2062 s , 0/ a i \u2062 s ) 1 / ( 1 + \gamma ),$

(A24)

$ c 1= ( a C I o \u2062 s , 0 ) n C 1 + ( a C I o \u2062 s , 0 ) n C,$

(A25)

Acknowledgments

We wish to thank Trevor Lamb for inspiring discussions and helpful comments on the manuscript. This research is supported by the Netherlands Organization for Scientific Research (NWO/ALW).

Commercial relationships: none.

Corresponding author: Hans van Hateren.

Email: j.h.van.hateren@rug.nl.

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

References

Brown, K. S.
(2000).

*Lead-lag algorithms.*. Retrieved from.
Dacey, D.
Packer, O. S.
Diller, L.
Brainard, D.
Peterson, B.
Lee, B.
(2000). Center surround receptive field structure of cone bipolar cells in primate retina.

*Vision Research*, 40, 1801–1811. [PubMed] [CrossRef] [PubMed]
Detwiler, P. B.
Hodgkin, A. L.
McNaughton, P. A.
(1980). Temporal and spatial characteristics of the voltage response of rods in the retina of the snapping turtle.

*Journal of Physiology*, 300, 213–250. [PubMed] [Article] [CrossRef] [PubMed]
Friedburg, C.
Allen, C. P.
Mason, P. J.
Lamb, T. D.
(2004). Contribution of cone photoreceptors and post-receptoral mechanisms to the human photopic electroretinogram.

*Journal of Physiology*, 556, 819–834. [PubMed] [Article] [CrossRef] [PubMed]
Hamer, R. D.
(2000). Computational analysis of vertebrate phototransduction: Combined quantitative and qualitative modeling of dark- and light-adapted responses in amphibian rods.

*Visual Neuroscience*, 17, 679–699. [PubMed] [Article] [CrossRef] [PubMed]
Hamer, R. D.
Nicholas, S. C.
Tranchina, D.
Lamb, T. D.
Jarvinen, J. L.
(2005). Toward a unified model of vertebrate rod phototransduction.

*Visual Neuroscience*, 22, 417–436. [PubMed] [Article] [CrossRef] [PubMed]
Hamer, R. D.
Nicholas, S. C.
Tranchina, D.
Liebman, P. A.
Lamb, T. D.
(2003). Multiple steps of phosphorylation of activated rhodopsin can account for the reproducibility of vertebrate rod single-photon responses.

*Journal of General Physiology*, 122, 419–444. [PubMed] [Article] [CrossRef] [PubMed]
Hood, D. C.
Finkelstein, m. A.
Boff, K. R.
Kaufman, L.
Thomas, J. P.
Sensitivity to light.

*Handbook of perception and human performance*. (pp. 5‐1––5‐66). New York: John Wiley.
Kelly, D. H.
(1961). Visual responses to time-dependent stimuli: 1 Amplitude sensitivity measurements.

*Journal of the Optical Society of America*, 51, 422–429. [CrossRef] [PubMed]
Kenkre, J. S.
Moran, N. A.
Lamb, T. D.
Mahroo, O. A. R.
(2005). Extremely rapid recovery of human cone circulating current at the extinction of bleaching exposures.

*Journal of Physiology*, 567, 95–112. [PubMed] [Article] [CrossRef] [PubMed]
Kennedy, M. J.
Dunn, F. A.
Hurley, J. B.
(2004). Visual pigment phosphorylation but not transducin translocation can contribute to light adaptation in zebrafish cones.

*Neuron*, 41, 915–928. [PubMed]. [CrossRef] [PubMed]
Lamb, T. D.
Pugh, Jr., E. N.
(2004). Dark adaptation and the retinoid cycle of vision.

*Progress in Retinal and Eye Research*, 23, 307–380. [PubMed] [CrossRef] [PubMed]
Lee, B. B.
Dacey, D. M.
Smith, V. C.
Pokorny, J.
(2003). Dynamics of sensitivity regulation in primate outer retina: The horizontal cell network.

*Journal of Vision*, 3, (7):5, 513–526, http://journalofvision.org/3/7/5/, doi:10.1167/3.7.5. [PubMed] [Article] [CrossRef]
Lee, B. B.
Pokorny, J.
Smith, V. C.
Martin, P. R.
Valberg, A.
(1990). Luminance and chromatic modulation sensitivity of macaque ganglion cells and human observers.

*Journal of the Optical Society of America*, 7, 2223–2236. [PubMed] [CrossRef] [PubMed]
Mahroo, O. A.
Lamb, T. D.
(2004). Recovery of the human photopic electroretinogram after bleaching exposures: Estimation of pigment regeneration kinetics.

*Journal of Physiology*, 554, 417–437. [PubMed] [Article] [CrossRef] [PubMed]
Nikonov, S.
Lamb, T. D.
Pugh, Jr., E. N.
(2000). The role of steady phosphodiesterase activity in the kinetics and sensitivity of the light-adapted salamander rod photoresponse.

*Journal of General Physiology*, 116, 795–824. [PubMed] [Article] [CrossRef] [PubMed]
Normann, R. A.
Perlman, I.
(1979). The effects of background illumination on the photoresponses of red and green cones.

*Journal of Physiology*, 286, 491–507. [PubMed] [Article] [CrossRef] [PubMed]
Press, W. H.
Teukolsky, S. A.
Vetterling, W. T.
Flannery, B. P.
(1992).

*Numerical recipes in Fortran*. New York: Cambridge University Press.
Pugh, E. N.
Lamb, T. D.
Stavenga,, D. G.
de, W. J.
Pugh,, E. N.
(2000).
Phototransduction in vertebrate rods and cones: Molecular mechanisms of amplification, recovery and light adaptation.

*Handbook of biological physics*. (3, pp. 183–254). Amsterdam: Elsevier.
Rushton, W. A.
Henry, G. H.
(1968). Bleaching and regeneration of cone pigments in man.

*Vision Research*, 8, 617–631. [PubMed] [CrossRef] [PubMed]
Stockman, A.
Langendrfer, M.
Smithson, H. E.
Sharpe, L. T.
(2006). Human cone light adaptation: From behavioral measurements to molecular mechanisms.

*Journal of Vision*, 6, (11):5, 1194–1213, http://journalofvision.org/6/11/5/, doi:10.1167/6.11.5. [PubMed] [Article] [CrossRef]
Tranchina, D.
Sneyd, J.
Cadenas, I. D.
(1991). Light adaptation in turtle cones.

*Biophysical Journal*, 60, 217–237. [PubMed] [Article] [CrossRef] [PubMed]
Valeton, J. M.
van Norren, D.
(1983). Light adaptation of primate cones: An analysis based on extracellular data.

*Vision Research*, 23, 1539–1547. [PubMed] [CrossRef] [PubMed]
van Hateren, J. H.
(2005). A cellular and molecular model of response kinetics and adaptation in primate cones and horizontal cells.

*Journal of Vision*, 5, (4):5, 331–347, http://journalofvision.org/5/4/5/, doi:10.1167/5.4.5. [PubMed] [Article] [CrossRef]
van Hateren, J. H.
Lamb, T. D.
(2006). The photocurrent response of human cones is fast and monophasic.

*BMC Neuroscience*, 7, 34. [CrossRef] [PubMed]
van Hateren, J. H.
Snippe, H. P.
(2006). Phototransduction in primate cones and blowfly photoreceptors: Different mechanisms, different algorithms, similar response.

*Journal of Comparative Physiology: A, Neuroethology, Sensory, Neural, and Behavioral Physiology*, 192, 187–197. [PubMed] [CrossRef] [PubMed]