We simulate a spatio-spectral photoreceptor mosaic, C, of size N_x × N_y × N_λ by assigning a cone spectral sensitivity of length N_λ in each location of an N_x × N_y spatial mosaic. Let (λmax,jM,λmax,jL) be the jth pair of peak sensitivity values drawn from a database of N_s pairs. Let A = (_{axy)Nx × Ny} be a matrix describing the spatial arrangement of M and L cones, where the elements are cone type labels: a_xy ∈ {M, L}. Let ψ(λ,λmax,jaxy) be a function that returns the value at wavelength λ of a spectral sensitivity with peak sensitivity at λmax,jaxy. To build the spatio-spectral mosaic for the jth spectral sensitivity pair, we extend the formulation in Sumner and Mollon (2000b) to include spatial dimensions as:

where x and y are the cone (pixel) coordinates; λ is wavelength sampling positions drawn from the N_λ-length vector of all wavelength sampling positions, Λ; and Δ(λ) gives the combined wavelength-dependent optical densities of the lens and macular pigment at wavelength λ. Thus, Equation (1) gives the spectral sensitivity at wavelength λ of a cone of type a_xy (either M or L), adjusted to account for the ocular media.

To simulate viewing of the scene, we present the system with a spectral image, P_i, of the same spatial and spectral dimensions as C_j (N_x × N_y × N_λ) and follow previous works (Kelber et al., 2003) by calculating quantal catch at each spatial location to produce an N_x × N_y quantal catch matrix, Q_ij, where:

It is assumed that the illuminant is already included with the spectral image. Notice that Equation (2) does not discriminate between cone types, so matrix Q_ij contains the quantal catches of both M and L cones.

We model the opponent processing of retinal ganglion cells (RGCs) as Difference of Gaussians (Soodak, 1986). Let W^ON be an ON-centre OFF-surround receptive field and W^OFF be an OFF-centre ON-surround receptive field. Each of these is formed as a two-dimensional circularly symmetrical Difference of Gaussians. Because of the symmetry, we define a 2D Gaussian, N(μ, σ), using a scalar mean, μ, and scalar standard deviation, σ, as:

The two receptive fields can then be defined as sums of zero-mean Gaussians with different standard deviations, as:

where σ_c and σ_s are the scalar standard deviations of the centre and surround Gaussians, respectively, and ω allows the sensitivity of the surround receptive field to be adjusted with respect to that of the centre.

Note that the receptive field does not discriminate between cone types, so this is a non-selective receptive field model. In order to achieve colour opponency, the centre must have high cone purity (draw input predominantly from one cone type). In the present work, we only model the single-cone centres of RGCs in the foveola, so cone purity of centres is 100%.

Next, we model the opponent processing of the ON-centre and OFF-centre RGCs. It is likely that midget RGCs cannot discriminate between L and M cones (Paulus and Kroger-Paulus, 1983; Benson et al., 2014), so let all ON-centre RGC responses be held in a N_x × N_y matrix, GijON. Moreover, let GijOFF be an equivalent matrix of OFF-centre RGC responses. These are computed as:

where ^* represents the convolution operation. Notice that both cone types (M and L) in Q_ij are involved in the convolution; i.e., this is a model of non-selective RGC centres and surrounds, so colour opponency is achieved by the fact that the centre contains only a single cone—the spectral sensitivity of the single-cone centre is exactly the same as that of its single cone, while the spectral sensitivity of the surround is a mixture of the spectral sensitivities of M and L cones.

The convolution operation in Equation (5) produces N_w − 1 pixels that incur border effects in GijON and GijOFF, where N_w is the length of each dimension of the window containing the receptive field. These are artefacts, they are removed prior to measuring performance.

The quantal catches of both cone types are contained within the same matrix, Q_ij, but subsequent steps treat M and L cones differently. From the given photoreceptor mosaic, we create two logical matrices, L and M, which act as pointers to cones of their respective type in Q_ij, as:

We now have two matrices of RGC responses: GijON contains both L ON-centre and M ON-centre (we will refer to these as LijON and MijON) responses, and GijOFF contains both L OFF-centre and M OFF-centre (LijOFF and MijOFF) responses. To build a red-green signal, we require: (i) only the RGCs that are excited by “reddish” stimuli; (ii) RGCs with both M and L centre so as to achieve full coverage of the scene. These two conditions are met by LijON and MijOFF. Assuming the cortex can learn RGC centre identities from the statistics of their activations (Wachtler et al., 2007), then LijON=LGijON and MijOFF=MGijOFF. So, the red-green matrix of outputs to the cortex, O_ij, is calculated as:

Similarly, the luminance signal, Oj′, is a combination of L and M centre RGC signals for full scene coverage, but only those with ON-centre receptive fields (which is equivalent to GijON). This is calculated as:

Note that the matrices O_ij and Oij′ are inputs to the cortex; therefore, no attempt has been made to model the spatial integration that is thought to occur in the cortex, but for which the mechanism has not yet been revealed (Solomon and Lennie, 2007).

Ideally, we would use the parameter settings (e.g., focal length, transmission of optical media) of trichromatic primates at the time when they diverged from dichromats. However, such information is not available. Unless otherwise stated, we use humans as models since the data is more comprehensive than for other species, and, where comparisons have been made, values tend to be similar across primate species (Cooper and Robson, 1969; Snodderly et al., 1984; Tovée et al., 1992).

For computational convenience, we model the photoreceptor mosaic as a square grid (Benson et al., 2014). This is a simplification of the hexagonal grid found in the primate eye (Hofer et al., 2005) that eases the computation and description, but does not affect our results because within one simulation we compare all spectral sensitivity pairs on exactly the same mosaic. For convenience and efficiency, we model one cone as equal in size to one pixel. In the human fovea, the minimal centre-to-centre spacing of cones is 1120 degrees visual arc, or 2.5 μm (Wandell, 1995), so we define one pixel as 1120×1120 degrees visual arc.

To demonstrate the advantage of spatio-chromatic modelling to estimate optimal spectral sensitivity tuning, we focus upon the spectral tuning of the cones that underlie the red-green system. For this reason, we only model the central region of the primate foveola, which contains no S cones. In humans, the S-cone-free region is around 20/60 degrees across (100 μm, or 40 cones) (Bumsted and Hendrickson, 1999), though it is smaller in macaques, at around 9/60 degrees across (45 μm, or 18 cones) (de Monasterio et al., 1985). However, varying the size of the mosaic does not affect any calculations, so we do not strictly limit ourselves to the size found in nature. Instead we choose the horizontal and vertical photoreceptor mosaic dimensions, N_x and N_y, to match the size of the spectral images.

For the L:M ratio of the photoreceptor mosaic, we use 1:1. This ratio has been observed as average in most non-human primates (Mollon and Bowmaker, 1992), so we assume the earliest trichromatic ancestor was similar. Though large variation has been observed in this ratio (Deeb et al., 2000), we do not include this in our simulation because it is not clear whether this variation existed in the earliest primate that evolved trichromacy and it is also not yet clear what mechanisms assure similar colour vision performance across different L:M ratios (Brainard et al., 2000).

To assign cone identities to the mosaic, let A = (_{axy)Nx × Ny} be the photoreceptor mosaic, and a_xy ∈ {M, L} be the cone identity at pixel (x, y). In the primate retina, L and M cones have a locally random arrangement (Mollon and Bowmaker, 1992). Thus, we randomly assign each a_xy with a 0.5 chance to be M, otherwise it is L.

For the combined wavelength-dependent optical densities of the lens and macular pigment, Δ, we follow Sumner and Mollon (2000a) in using those of humans, as given in Wyszecki and Stiles (1982). For the function which generates the spectral sensitivities from the peak sensitivities, ψ(·, ·), we use the peak wavelength-dependent nomogram described in Stockman and Sharpe (2000) and the resulting absorbance curves are converted to absorptance curves.

The size of the receptive fields is chosen to match that in the primate fovea in having a single cell centre (Calkins and Sterling, 1996), which is achieved by σ_c = 0.25 cones (95% of the volume of a Gaussian falls within 2 standard deviations either side of the mean, and 4σ_c = 1). The size of the surround is determined from the centre-to-surround ratio as: σ_c/σ_s = 0.15 (Croner and Kaplan, 1995), which yields σ_s = 1.66 cones (95% of the Gaussian is covered by a diameter of 4σ_s = 6.67 cones). N_w, the length of each dimension of the window that contains the receptive field, is set to 9 cones. This accommodates slightly less than three standard deviations each side of the mean (6σ_s = 9.96), but allows us to keep a single cone in the centre of the window. For our fixed receptive field size, we verified that this window size did not influence results by running a comparison test with N_w = 11 cones and observed no difference in the results. It has been shown that the relative sensitivities of centre and surround are highly variable from cell to cell, but that the sensitivity of the surround is, on average, 55% weaker than that of the centre (Croner and Kaplan, 1995). We use this average value by setting ω = 0.55.

We use grating patterns rather than images of natural scenes because this allows us to vary spatial frequency in a well-defined manner. To emulate the sinusoidal gratings commonly used in psychophysical studies (Mullen, 1985; Brainard et al., 2008; Lee et al., 2012), we create two-colour grating patterns containing ecologically valid reflectance spectra in each pixel (spectral images). Each spectral sensitivity pair's performance is measured as its ability to maximise the distinctiveness of the spatial pattern; i.e., the sinusoidal grating. The highest spatial frequency used is 4 cycles per degree, as it has been shown that significant chromatic aberration affects the retinal image at spatial frequencies above 4 cycles per degree (Flitcroft, 1989). The exact spatial frequencies, in cycles per degree, are {4, 2, 1, 0.5}, which equates to {30, 60, 110, 222} cones (1 cone = 1 pixel) per cycle, or {1560,3060,5560,11160} degrees visual arc, respectively. The creation of a spectral image is demonstrated in Figure 3, and other spectral images vary in the spatial frequency of the grating (A), the fruit spectrum (B), and the set of luminance coefficients (D).

All spectral images in one simulation use the same spatial pattern: P¯ = (p¯xy)Nx×Ny,p¯xy∈[0, 1], which is also used as the reference pattern with which the model outputs are compared. For each simulation, we create the database of spectral images, Pi∈ℝNx×Ny×Nλ, i=1,2,…,Nf. Let si∈ℝNλ be the ith target spectrum in the database of N_f target spectra, and let sb∈ℝNλ be the mean of all background leaf spectra. For each pixel position (x, y) in P_i, the spectrum is composed from the target and background incident spectra, s_i and s_b, as  p¯xysi+(1- p¯xy)sb.

The spectra are taken from the Cambridge database of natural spectra (Regan et al., 1998; Sumner and Mollon, 2000b), which contains a variety of spectra recorded from fruits eaten by several primate species, and mature leaves from the same environments. The spectra are provided at 4 nm sampling intervals, and we truncate them to the primate visual range of 400–700 nm, so N_λ = 76 and Λ = {400, 404, 408, …, 700}. The number of images in the database matches the number of spectra in the Cambridge database: N_f = 1139 fruits. Note that each spectral image, P_i, is created from only two spectra: one target and one background, where the background spectrum is always the mean of all (409) mature leaves, as has been done in previous work (Chittka and Menzel, 1992; Osorio and Vorobyev, 1996). We also tested using randomly selected background leaf spectra in each pixel rather than the mean, but this had no effect on the results. For one set of simulations, we use spectra of Munsell chips (Parkkinen et al., 1989). In this case, there is no constant background spectrum for all images. Instead, for each image, we draw two spectra from the dataset randomly and without replacement.

Prior to building the images, the illuminant was added to all spectra as: s = i ⊙ s′, where s′ is the original spectrum, i∈ℝNλ is the illuminant, and ⊙ is the element-wise multiplication operator. We tested using three different illuminants: the D65 standard illuminant and two for forest areas in French Guiana and Uganda. We found that all three illuminants resulted in the same set of optimal spectral sensitivities, and so have selected to present only the results produced using the Ugandan forest illuminant from the Cambridge dataset (Regan et al., 1998; Sumner and Mollon, 2000b).

In natural scenes, there is far more variance in luminance compared to chromaticity (Ruderman et al., 1998; Wachtler et al., 2001), mainly due to the effects of shadows. Such effects are characterised as multiplicative factors on the reflected spectrum (Rubin and Richards, 1982). Though natural scenes may also contain some additive effects due to reflectance, we assume that these would be small enough to be negligible, as the specularities they cause would only occur when the surface was at a certain angle relative to the eye, and the animal could move its head to change this. Therefore, we do not include any additive reflectance effects. We included luminance variation by creating greyscale images which capture the 1/f power spectrum typically observed in natural images (Field, 1987; Billock, 2000; Millane et al., 2003), where f is the spatial frequency, and then using the pixel values as coefficients for the spectra in our images (1/f noise is also known as Pink noise). We used a set of 49 natural leafy scenes gathered from an internet search with the constraints that the images must contain only leaves, branches and sky, with branches and sky being minimal. These images were all converted to greyscale. We then performed a Fast Fourier Transform (FFT) on each image and used the median in each pixel value to create average magnitude and phase images. We verified that the resulting power spectrum was a close fit with the 1/f slope, then these median values were transformed back into the spatial domain via inverse FFT. An example luminance coefficient image is shown in Figure 3D.

The wavelength-dependent refractive effect of the lens is, strictly, a part of the biological model rather than the image model. However, as this is constant across all images and regardless of the cone spectral sensitivities, we include it in the spectral images. To account for the refractive effect of the lens, we apply a pixelwise point-spread function to all spectral images by convolving each wavelength layer of the spectral images with a wavelength appropriate Airy disk. The Airy disk is created using the data available for humans: lens refractive index = 1.406 (Garner et al., 1998), focal length = 21.3 mm (Van Norren and Tiemeijer, 1986), pupil diameter = 5 mm (Liang and Williams, 1997), and for each wavelength λ ∈ Λ. The line spread functions for the shortest and longest wavelengths we use, λ = 400 nm and λ = 700 nm, are shown in Figure 3F.

Vertical image size is N_y = 30 cones (pixels). Horizontal image size, N_x, is dependent on spatial frequency, such that each image contains two cycles of the sinusoidal grating. Specifically, this gives N_x ∈ {60, 120, 220, 444} cones (pixels).

We measure the ability of each peak sensitivity pair (λmax,jM,λmax,jL), to facilitate the discrimination of chromatic gratings. We make the assumption that the cortical cells receiving input from the parvocellular channel attempt to estimate the input image (Manning and Brainard, 2010). The jth peak sensitivity pair's estimate of the ith image is stored in matrices of outputs to the cortex: O_ij for the red-green signal and Oij′ for the luminance signal. Moreover, as we simulate the spectral images, the reference pattern, P¯, is known, so the performance is measured as the similarity between the reference pattern and the matrices of outputs to the cortex. Here, we will describe how the performance measure is calculated for the red-green signal by comparison of O_ij and P¯. Performance of the luminance signal is measured in the same manner—simply replace O_ij with Oij′ and z_j with zj′.

After finding O_ij for all i = 1, 2, …, N_f images, we compare them all with the reference pattern, P¯, and produce a single measure of the performance, z_j, of spectral sensitivity pair j. Values in P¯ vary from 1 in target locations to 0 in background locations, while O_ij contains positive values in cells which are redder than their local neighbourhood (determined by the size of the receptive field surround) and negative values in cells which are greener than their local neighbourhood. To ensure O_ij and P¯ are in the same range, we normalise O_ij to [0, 1] prior to the comparison. As all spectral images in one simulation have the same spatial pattern, P¯ is the same for all P_i. To measure the similarity between a single O_ij and P¯, we use the Peak Signal-to-Noise Ratio, which is commonly used in image processing for comparing a reference and filtered image (Drew and Bergner, 2007). Its calculation is based upon Mean Square Error, but taking into account the maximal possible signal power. As such, this performance is maximal for the spectral sensitivity pair that causes the greatest mean distance in red-green space between all fruits and the background leaves. The performance measure, z_j, is calculated as:

where PSNR(·, ·) is the is the Peak-Signal-to-Noise Ratio of the two arguments, defined as:

where MSE(Oij,P¯) is the Mean Square Error of its two arguments. In Equation (10), peakval = 1 is the peak value that the signal can take, and the logarithm converts the result to decibels. The spatial performance is implicit because the red-green value at a given location in O_ij is a function of the activations of all cones in a local neighbourhood (determined by the size of the RGC receptive field). O_ij will contain noise which arises because: (i) the receptive field surrounds synapse both M and L cones (they are non-selective), and the ratio of L:M cones differs from surround to surround due to the random mosaic arrangement; (ii) the luminance variation can cause the activations of the cones stimulated by the spectra of the same material to be less correlated than cones stimulated by the spectra of different materials, particularly for very similar spectra. Noise of the first type should be minimised by smaller λmaxL-λmaxM, while noise of the second type will favour some separation, with the amount being dependent upon the nature of the spectra. The best performing spectral sensitivity pair is the one that maximises the chromatic signal described by P¯ while minimising noise.

We aim to determine whether a spatio-chromatic parvocellular channel model will predict highly overlapping M and L spectral sensitivities even if a spectrally rich model of the environment is used. This is a control to ensure that the spatio-chromatic model is not biased toward predicting spectral sensitivities similar to those found in primates. We use a database of spectra of Munsell chips (Parkkinen et al., 1989) to create the spectral images, vary both spectral sensitivities over the range: λmax,jM,λmax,jL∈{490, 500, 510, …, 550, 560}, and apply the constraint λmax,jL≥λmax,jM. Each simulation is repeated 50 times on differently seeded random mosaics.

The results show that discrimination performance increases the more separated λmaxM and λmaxL become, as shown in the performance plot in Figure 4A. Accordingly, the optimal peak spectral sensitivities are: λmaxM= 490 nm and λmaxL= 560 nm. Performance plots at all spatial frequencies were qualitatively similar, so only the result for 4 cycles per degree gratings is shown.

We aim to determine whether a spatio-chromatic model requires a long-wavelength limit of λmaxL≤ 562 nm to constrain the L cone spectral sensitivity to its observed position, rather than a longer wavelength. We assume the task is foraging for fruit and perform an exhaustive search where both peak sensitivities, λmaxM and λmaxL, vary. We apply the constraint that λmaxL≥λmaxM. The minimum and maximum of 490 and 598 nm cover the range over which mammal LWS opsins exist (Jacobs, 2009), with an additional ~40 nm in the longwave direction to investigate whether the optimal L cone peak would lie at a longer wavelength than found in nature if no biological constraints applied, i.e., λmax,jM,λmax,jL∈{490, 494, 498, …, 594, 598}. Each simulation is repeated 50 times on differently seeded random mosaics.

The results are presented in Figure 4B. The performance plots at all spatial frequencies were qualitatively similar, so only the result for 4 cycles per degree gratings is shown. On all simulations, λmaxL occupies the longest wavelength possible, and for any given λmaxM, performance is always better if λmaxL> 562 nm. When both λ_max are at relatively short wavelengths (i.e., when the mean of the two is shorter than ~540 nm), discrimination performance is lower than the nearest configuration for which λmaxM=λmaxL (at which point this ceases to be a red-green signal and becomes a luminance signal). The absolute optimum configuration occurs at λmaxL= 598 nm and λmaxM= 521.6 nm, whereas when the long-wavelength limit is applied, the optimum λmaxM is 525.8 nm. The extreme optimal λmaxM on individual runs were 518 and 526 nm—this variation is due to the random “seeding” of the photoreceptor mosaic between simulation repetitions and, to a lesser degree, of the luminance coefficients between images. These results guided us to limit the search range for subsequent simulations sets to 515 ≤λmaxM≤ 535 nm in 1 nm increments and λmaxL= 562 nm (discrimination performance is always better for λmaxL at longer wavelengths, so this is equivalent to λmaxL≤ 562 nm).

We aim to determine whether the optimal position of the M cone spectral sensitivity is more accurately predicted by a spatio-chromatic analysis than by a purely spectral analysis in which spatial dimensions are not considered. We assume the ecological task of foraging for fruit, as this facilitates comparison with the purely spectral analyses of previous works (Osorio and Vorobyev, 1996; Regan et al., 1998; Sumner and Mollon, 2000a; Lewis and Zhaoping, 2006). We use the results of the long-wavelength limit simulations to limit our search space, thereby facilitating a higher-resolution (smaller increment between λ_max positions) search. As λmaxL always occupies the longest wavelength available, fixing λmaxL= 562 nm is equivalent to λmaxL≤ 562 nm. The search interval is reduced to 1 nm, and we increase the number of repetitions with differently seeded mosaics to 100. As our model produces both red-green and luminance signals that are sent to the cortex, we also demonstrate the spatio-chromatic advantage of the red-green system over the luminance system for the task of fruits foraging.

The results are presented numerically in the upper row of Table 2 and visually in Figure 5. For the red-green system (Figure 5A), the optimal λmaxM= ~525 nm, with standard deviations of ~1 nm and a 95% confidence interval of <0.5 nm at all spatial frequencies. There is a small but significant change in λmaxM with spatial frequency [p = 8.9 × 10⁻¹², α = 0.01, two-sample t-test between data for highest (λmaxM= 525.8 nm) and lowest (λmaxM= 524.84 nm) spatial frequencies]. The red-green system performs relatively better than the luminance system (Figure 5B) at all λmaxM positions in the range we tested. For both systems, discrimination performance increases as spatial frequency decreases. Means, standard deviations and confidence intervals are not given for the luminance system as the optimal λmaxM always occurs at the longest wavelength available. A note on the results of the luminance system: while it is typical to observe size-tuned responses in the luminance system, our results do not show this. Instead, the tuning curves for lower spatial frequencies have higher z scores. It has been observed that contrast sensitivity and the band-pass shape decrease as luminance levels decrease (Enroth-Cugell and Robson, 1966). As our goal here is to provide a comparison for the red-green system, we have used the same spectral images, and these contain fruit and leaf spectra are highly similar in intensity. Thus, the higher performance at lower spatial frequencies is a result of the low contrast of red-green gratings to the luminance system.

We aim to determine the effect that luminance variation has on the position of the optimal λmaxM. luminance variation is a source of spatial noise—i.e., it causes different regions of the cone image to be active to different degrees even when stimulated by the same material. As such, it is a noise source that could not be included in previous works that did not model spatial dimensions. We applied the long-wavelength limit by fixing λmaxL= 562 nm, and used a broader but coarser search range for λmaxM of 440–560 nm as during testing we observed a shift to shorter wavelengths. As our aim with this test is only to identify qualitative shifts in the optimal λmaxM due to luminance variation, we reduce the sampling interval to 10 nm. To remove the luminance variation, the luminance coefficients are all set to 1 when creating the spectral images.

The results are presented numerically in the bottom row of Table 2 and visually in Figure 6. Without luminance variation, λmaxM in the range 440–500 nm lead to similar performance, but at longer wavelengths than ~500 nm, performance decreases with wavelength. There is a significant effect of spatial frequency, with the optimal λmaxM occurring at a longer wavelength for higher spatial frequencies (p = 5.2 × 10⁻³¹, α = 0.01, two-sample t-test between data for highest and lowest spatial frequencies: λmaxM= 483.68 nm and λmaxM= 444.38 nm, respectively).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.