Simulating human cones from mid-mesopic up to high-photopic luminances

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 ⁻³ to 10 cd/m ²) to photopic (10 to >10 ⁶ cd/m ²), 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⁶ 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

 

(A) Model of the human cone. Boxes containing a τ represent first-order low-pass filters of unit DC gain; other boxes represent linear or nonlinear transformations of variables. Units and notation are consistent with van Hateren (2005). The gain for photon capture is proportional to the fraction of excitable pigment (1 − B − c_NR*, where B is the bleached pigment and R* is the excited pigment). R* has a lifetime τ_R before it is bleached to B which has an effective lifetime τ_B before it is reconverted to excitable pigment. R* is multiplied by a scaling constant c_N to obtain normalized units consistent with Mahroo and Lamb (2004). R* excites a G protein that quickly forms a complex E* (with lifetime τ_E) with phosphodiesterase (PDE). PDE can hydrolyze cGMP at a rate β which consists of a dark activity c_β by unexcited PDE and a light-dependent part k_βE*. Saturation of cGMP hydrolysis reduces β to an effective β_e according to a compressive nonlinearity. The nonlinear differential equation describing cGMP hydrolysis produces a gain factor 1/β_e and a time constant τ_X = 1/β_e, implicitly depending on light intensity (van Hateren, 2005). The cGMP concentration (X) opens transduction channels (with apparent cooperativity n_X), producing current I_os into the cone's outer segment. Part of the current consists of calcium (C extruded with a time constant τ_C), which drives cGMP synthesis (with cooperativity n_C) and, thus, produces a feedback gain α. I_os drives the nonlinear membrane of the inner segment, with a gain dependent on the membrane potential V_is. The model neglects a small overall delay of 1–3 ms (van Hateren, 2005; van Hateren & Lamb, 2006), which may be partly due to diffusional processes not explicitly modeled. (B) PDE saturation as an alternative to saturating cGMP hydrolysis. NL: nonlinearity. (C) Regulation of pigment phosphorylation by calcium as an alternative to saturating cGMP hydrolysis. DNL: dynamic nonlinearity.

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_NR*, 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 I produces excited pigment R* according to the following differential equation:  

${\tau }_R\dot{R}*=I(1-B-c_{\mathrm{N}}R*)-R*,$

(1)

where the dot denotes time differentiation. Note that 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): 

$\dot{B}=c_{\mathrm{N}}R*/{\tau }_R-\frac{K_B}{{\tau }_{B,0}}\frac{B}{B+K_B}.$

(2)

The first term on the right-hand side represents the production of 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 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  

$\begin{array}{cc}& {\tau }_B\dot{B}=g_{B}R*-B,\hfill \\ \\ & \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{\tau }_B={\tau }_{B,0}\frac{B+K_B}{K_B}\mathrm{a}\mathrm{n}\mathrm{d}{g}_B=c_{\mathrm{N}}{\tau }_B/{\tau }_R.\hfill \end{array}$

(3)

This equation shows that the slow gain control provided by 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  

$0=I_0(1-B_0-c_{\mathrm{N}}{R}_0*)-R_0*$

(4)

or  

$R_0*=(1-B_0)\frac{I_0}{1+c_{\mathrm{N}}{I}_0}.$

(5)

The steady state of Equation 2 produces a quadratic equation in B ₀ when the R ₀* 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 ₀ ≪ K _B τ _R/ τ _B,0 or I ₀ ≪ 6,600 td. Then, it can be shown that R ₀ ≈ 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 ₀ ≫ K _B τ _R/ τ _B,0 or I ₀ ≫ 6,600 td, where Equation A22 reduces to  

$B_0=1-\frac{{\tau }_R}{{\tau }_{B,0}}\frac{K_B}{1+K_B}\frac{c_{\mathrm{N}}{I}_0+1}{c_{\mathrm{N}}{I}_0},$

(6)

and consequently (using Equation 5, R _N,0* = c _N R ₀*, and R _N,0 = 1 − B ₀ − R _N,0*)  

$R_{\mathrm{N},0}*=\frac{{\tau }_R}{{\tau }_{B,0}}\frac{K_B}{1+K_B}$

(7)

and  

$R_{\mathrm{N},0}=\frac{{\tau }_R}{{\tau }_{B,0}}\frac{K_B}{1+K_B}\frac{1}{c_{\mathrm{N}}{I}_0}.$

(8)

This implies that the excited pigment R _N,0* goes to a fixed fraction of about 2.3 × 10 ⁻⁵ (∼1,100 molecules) at high intensities, whereas the available excitable pigment R _N,0 becomes inversely proportional to I ₀. 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 ₀ ≫ 10 ⁶ 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 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 }_{\mathrm{e}}=\frac{\beta }{1+\beta /{\beta }_{\mathrm{e}\mathrm{,}\max }},$

(9)

where β _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 τ

$\dot{y}$

= 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\left(n\right)=f_{1}y(n-1)+f_{2}x(n-1)+f_{3}x\left(n\right),$

(10)

with 

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

(11)

where τ 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

 

Dark gray line: data from Figure 5C of Lee et al. (2003), showing the response (in millivolts per troland, average of 15 H1 horizontal cells) to a 19.5-Hz sinusoid of low contrast superimposed on various background intensities (abscissa). Dashed blue line: extrapolation following Weber's law (slope: −1). Solid red line: fit to the data and extrapolation with the H1 model of van Hateren (2005) extended with pigment bleaching and saturating cGMP hydrolysis. Fitted parameters: k_β = 1.4 × 10⁻⁴, β_e,max = 4, a_C = 0.23, a_is = 2.9 × 10⁻² (see the 1 for units and default values of the other parameters). Dashed red line: model predictions at 0.61 Hz.

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 ⁴ td, the curve obtains a slope of −1 and, thus, conforms to Weber's law (see also Equations 8 and A13). 

Figure 3

 

Responses to a stimulus as in Figure 2 at various stages in the model of Figure 1, normalized at log-td = 0. Pigment bleaching dominates at intensities >10 ⁴ td (curve marked R*); saturation of cGMP hydrolysis contributes between 10 ³ and 10 ⁴ td ( β _e), and below 10 ³ td hydrolysis of cGMP (1/ β _e) combines with calcium feedback ( α) to produce the current I _os and membrane potential V _is.

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 ³- to 10 ⁴-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 ³- to 10 ⁴-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 ²⁺ 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² or more, for background intensities of 1, 30, 100, 300, 1,000, 3,000, 10⁴, 3 × 10⁴, and 10⁵ td. The background intensity of each black curve coincides with the point where it crosses the red curve. 

Figure 4

 

Contrast responses at different light levels. The inset shows a response to a 100-ms flash to 200 td from a background of 100 td; normalization is such that the response is 0 in the dark and −1 for I_os = 0. The peak response to similar flashes is shown inverted in the main panel. The red trace denotes the response to a steady background, and the black traces indicate the digression from the steady state when negative or positive contrast flashes are given, for a range of flash intensities (abscissa). Black traces are for background intensities of 1, 30, 100, 300, 1,000, 3,000, 10⁴, 3 × 10⁴, and 10⁵ td.

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⁴ 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

 

Small-signal frequency response of cones, depending on background intensity (numbers at curves, in trolands). At high intensities and not too low frequencies (gray area), sensitivity is approximately inversely proportional to intensity (Weber's law); the high-frequency cutoff is changing considerably in the 1- to 1,000-td range; the low-frequency falloff at high intensities is a consequence of bleaching kinetics.

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 ⁵- and 10 ⁶-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 ⁶ 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⁶-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): 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 ²⁺(−); 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).  

$\frac{\mathrm{d}R*}{\mathrm{d}t}=\left[I\right(1-B-c_{\mathrm{N}}R*)-R*]/{\tau }_R$

(A1)

 

$\frac{\mathrm{d}B}{\mathrm{d}t}=c_{\mathrm{N}}R*/{\tau }_R-\frac{K_B}{{\tau }_{B,0}}\frac{B}{B+K_B}$

(A2)

 

$\frac{\mathrm{d}E*}{\mathrm{d}t}=[R*-E*]/{\tau }_E$

(A3)

 

$\frac{\mathrm{d}X}{\mathrm{d}t}=\alpha -{\beta }_{\mathrm{e}}X$

(A4)

 

$\frac{\mathrm{d}C}{\mathrm{d}t}=(I_{\mathrm{o}\mathrm{s}}-C)/{\tau }_C$

(A5)

 

$\frac{\mathrm{d}{V}_{\mathrm{i}\mathrm{s}}}{\mathrm{d}t}=(I_{\mathrm{o}\mathrm{s}}/g_{\mathrm{i}}-V_{\mathrm{i}\mathrm{s}})/{\tau }_{\mathrm{m}}$

(A6)

 

$\frac{\mathrm{d}{g}_{\mathrm{i}}}{\mathrm{d}t}=(a_{\mathrm{i}\mathrm{s}}{V}_{\mathrm{i}\mathrm{s}}^{\gamma }-g_{\mathrm{i}})/{\tau }_{\mathrm{i}\mathrm{s}}$

(A7)

where  

$\beta =c_{\beta }+k_{\beta }E*$

(A8)

 

${\beta }_{\mathrm{e}}=\frac{\beta }{1+\beta /{\beta }_{\mathrm{e}\mathrm{,}\max }}$

(A9)

 

$\alpha =1/[1+(a_{C}C{)}^{n_C}]$

(A10)

 

$I_{\mathrm{o}\mathrm{s}}=X^{n_X}$

(A11)

The following parameters and default values are used: c _N = τ _R K _cone/ R _tot = 4.1 × 10 ⁻⁹ (td) ⁻¹, with K _cone = 60 isom s ⁻¹ (td) ⁻¹ and R _tot = 5 × 10 ⁷ 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⁻³ (ms)⁻¹ (dark PDE activity); k_β = 1.4 × 10⁻⁴ (ms)⁻¹ (td)⁻¹ (E* dependence of PDE activity); β_e,max = 4 (ms)⁻¹ (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²⁺ 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⁻² (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  

$\begin{array}{l}I=I_0+i\\ R*=R_0*+r*\\ B=B_0+b.\end{array}$

(A12)

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 

$\begin{array}{c}\frac{\tilde{r}*}{\tilde{i}}=\frac{(1-B_0-c_{\mathrm{N}}{R}_0*)}{1+c_{\mathrm{N}}{I}_0+\frac{c_{\mathrm{N}}{I}_0}{\frac{K_B^2{\tau }_R}{(B_0+K_B{)}^2{\tau }_{B,0}}+i\omega {\tau }_R}+i\omega {\tau }_R}\\ \approx \frac{R_{\mathrm{N},0}}{1+i\omega {\tau }_R},\hfill \end{array}$

(A13)

where the second line of the equation follows from c_NI₀ ≪ 1, with the exception of very low ω (see Figure 5). c_NI₀ ≪ 1 is true for normal intensities because 1/c_N > 10⁸ 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  

$\begin{array}{l}E*=E_0*+e*\\ {\beta }_{\mathrm{e}}={\beta }_{\mathrm{e},0}+\mathrm{\Delta }{\beta }_{\mathrm{e}},\end{array}$

(A14)

Equation A3 yields  

$\frac{\tilde{e}*}{\tilde{r}*}=\frac{1}{1+i\omega {\tau }_E},$

(A15)

and a small-signal analysis of Equations A8 and A9 gives  

$\frac{{\widetilde{\mathrm{\Delta }\beta }}_{\mathrm{e}}}{\tilde{e}*}=\frac{k_{\beta }{\beta }_{\mathrm{e}\mathrm{,}\mathrm{\backslash max}}^2}{({\beta }_{\mathrm{e}\mathrm{,}\mathrm{\backslash max}}+{\beta }_0{)}^2}.$

(A16)

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 

$H_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{e}}\left(\omega \right)=H_R\left(\omega \right)H_E\left(\omega \right)H_{\beta }\left(\omega \right)H_{\beta \mathrm{,}\mathrm{o}\mathrm{s}}\left(\omega \right)H_{\mathrm{i}\mathrm{s}}\left(\omega \right),$

(A17)

with 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 ,\mathrm{o}\mathrm{s}}=\frac{-n_X}{{\beta }_{\mathrm{e}\mathrm{,}0}+{\beta }_{\mathrm{e}\mathrm{,}0}{n}_X{n}_C{c}_1/(1+i\omega {\tau }_C)+i\omega }$

(A18)

 

$H_{\mathrm{i}\mathrm{s}}=\frac{(1+i\omega {\tau }_{\mathrm{i}\mathrm{s}})V_{\mathrm{i}\mathrm{s},0}}{(1+i\omega {\tau }_{\mathrm{i}\mathrm{s}})(1+i\omega {\tau }_{\mathrm{m}})+\gamma }.$

(A19)

In Equation A18, the β 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₀ = I_os,0): 

$I_{\mathrm{o}\mathrm{s},0}^{1/n_X}+a_C^{n_C}{I}_{\mathrm{o}\mathrm{s},0}^{n_C/n_X}=1/{\beta }_{\mathrm{e}\mathrm{,}0},$

(A20)

with 

$\begin{array}{l}{\beta }_{\mathrm{e}\mathrm{,}0}=\frac{{\beta }_0}{1+{\beta }_0/{\beta }_{\mathrm{e}\mathrm{,}\max }}\\ {\beta }_0=c_{\beta }+k_{\beta }{E}_0*\\ E_0*=R_0*\\ R_0*=(1-B_0)\frac{I_0}{1+c_{\mathrm{N}}{I}_0}.\end{array}$

(A21)

Here, I₀ is the steady-state intensity (in trolands). B₀ follows from eliminating R₀* from the steady state of Equations A1 and A2, which yields 

$\begin{array}{cc}{B}_0& =[1-K_B-\frac{{\tau }_R{K}_B(c_{\mathrm{N}}{I}_0+1)}{{\tau }_{B,0}{c}_{\mathrm{N}}{I}_0}\hfill \\ \\ & +\left(\right(1-K_B-\frac{{\tau }_R{K}_B(c_{\mathrm{N}}{I}_0+1)}{{\tau }_{B,0}{c}_{\mathrm{N}}{I}_0}{)}^2+4K_B{)}^{\frac{1}{2}}]/2.\hfill \end{array}$

(A22)

The R_N,0 needed in Equation A13 then follows from 

$R_{\mathrm{N},0}=1-B_0-c_{\mathrm{N}}{R}_0*.$

(A23)

Once I_os,0 is obtained from Equation A20, one finds 

$V_{\mathrm{i}\mathrm{s},0}=(I_{\mathrm{o}\mathrm{s},0}/a_{\mathrm{i}\mathrm{s}}{)}^{1/(1+\gamma )},$

(A24)

as needed in Equation A19, and 

$c_1=\frac{(a_C{I}_{\mathrm{o}\mathrm{s},0}{)}^{n_C}}{1+(a_C{I}_{\mathrm{o}\mathrm{s},0}{)}^{n_C}},$

(A25)

as needed in Equation A18. 

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: [e-mail]. 

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