Spatial summation of individual cones in human color vision

Introduction

A central goal of vision science is to understand how signals from photoreceptors are transformed into sight and the limitations each stage of processing imposes on perception. Photoreceptors provide an organism with real-time information about the environment. However, the signals conveyed by individual neurons are noisy and ambiguous. One strategy for reducing uncertainty is to pool signals across multiple detectors. Under low light conditions, for example, the visual system combines signals from many hundreds of rod and cone photoreceptors in order to boost sensitivity . One drawback of signal pooling is a loss in spatial resolution: both acuity and contrast sensitivity are reduced under low-light levels . We studied the influence of spatial pooling on the color appearance of cone-targeted spots.

Results

Load data and import libraries.

First load in data and set some parameters that will be used during plotting later.

# setwd("~/R/Schmidt-Boehm-Tuten-Roorda_2019")
library(tidyverse)
library(ggthemes)
library(gridExtra)
library(corrplot)
library(car)
# set the colors that will be used in plots below.
red <- "#e84040"
green <-"#40c631"
blue <- "#4042e8"
yellow <- "#e8df40"
gray <- "#595959"
# make sure a figures directory exists
dir.create('figures', showWarnings=FALSE)

trial_data <- read_csv('filtered_data.csv')
head(trial_data)

Example session

In a previous experiment, the one and two cone conditions were equated for detectability. We subsequently collected appearance measurements. All other stimulus conditions were identical to the detection task. Three cones were selected for study in each session (Fig~A). On each trial either a single cone or a pair was targeted. After each flash, the subject judged the color of the spot using a hue and saturation scaling paradigm . Each cone and pair was tested twelve times (72 trials per session). A total of 198 pairs were tested across three subjects. Hue and saturation scaling data were transformed into a color opponent representation. For each trial, the degree of perceived greenness versus redness and yellowness versus blueness was computed from percentage ratings as follows: \(gr = (green\% - red\%) / 100\%\) and \(yb = (yellow\% - blue\%) / 100\%\). In this representation, saturation is expressed as the distance from the origin (in city block metric). A pure white response falls at the origin and a maximally saturated report falls along the outer diamond.

The results of one session are plotted in Fig~B. In this example, Cone 1 was an M-cone and had a bias towards green. Cone 2 was an L-cone and elicited predominantly white reports. Cone 3, also an L-cone, was rated reddish-yellow (orange) with medium saturation. The motivation for this study was to understand the algorithm the visual system uses when combining information across cone pairs. In the example, when Cone 1 was targeted together with either Cone 2 or Cone 3, the average report was desaturated. In comparison, when Cone 2 and 3 were targeted they elicited a medium saturated orange report. Below, we analyze the results from all sessions and subjects.

example <- trial_data %>%
  filter(., subject == '20076R', session== 4) %>%
  group_by(., masterID1, masterID2, isPair) %>% 
  summarise(., yb_mean=signif(mean(yb), 3), 
            yb_sem=signif(sd(yb) / sqrt(length(yb)), 3), 
            gr_mean=signif(mean(gr), 3), 
            gr_sem=signif(sd(gr) / sqrt(length(gr)), 3), N=length(gr),
            typeIDstr=typeIDstr[1])
boundary <- data.frame(x=c(-1, 0, 1, 0, -1), y=c(0, 1, 0, -1, 0))
(examplePlot <- ggplot(data=example, 
                       aes(x=yb_mean, 
                           y=gr_mean,
                           color=interaction(typeIDstr, isPair),
                           shape=isPair)) +
    geom_point(size=3) + 
    geom_errorbar(aes(x=yb_mean, ymin=gr_mean - gr_sem, ymax=gr_mean + gr_sem)) +
    geom_errorbarh(aes(y=gr_mean, xmin=yb_mean - yb_sem, xmax=yb_mean + yb_sem, height=0.01)) + 
    coord_equal(xlim=c(-1, 1), ylim=c(-1, 1)) + 
    xlab('(yellow - blue) / total') + 
    ylab('(green - red) / total') + 
    scale_color_manual(
      values = c(red, green, red, yellow), 
      name="spectral type") +
    geom_text(label=c('1', '1+2', '1+3', '2', '2+3', '3'),
              nudge_x=-0.15) +
    geom_path(data = boundary, aes(x=x, y=y), colour='gray', linetype=2, inherit.aes = F) +
    theme_bw(base_size = 15) + 
    theme(legend.position="none")
  )

ggsave("figures/example_session.pdf")

Variability in sensations from cones with the same sensitivity

The motivating question behind these experiments was how does the visual system combine inputs from cones when making color appearance judgments? To begin to answer this question, we grouped each trial based on which cone or pair was probed. The results are reported in a table below. In both subjects with classified cones, individual and pairs of M-cones produced gr means that were greater than zero, while L-cone conditions led to negative values. This observation is consistent with a predictive relationship between cone type and color report, as found previously (Sabesan et al. 2006; Schmidt et al. 2018).

Each point in these plots represents the mean response measured from a single cone or a pair. This plot illustrates the variability in responses across cones/pairs and between subjects. There are a few features to note. Firstly, within a single subject, there was considerable variability between cones and pairs with the same spectral sensitivity. Similar variability in sensations from single cones has been reported previously (Hofer et al. 2005; Sabesan et al. 2016; Schmidt et al. 2018). This is the first report of variability in sensations elicited from pairs of cones. Secondly, there were individual differences in color responses: S20075 used blue more frequently than the two other subjects and S10001 did not report yellow on any trials. However, the general patterns are similar. Most of the variance was found along the green-red dimension; there were few points that fell in the blueish-red or greenish-yellow quadrants. In the two subjects with classified mosaics, we additionally found L-cone targeted trials tended produce a reddish bias, while M-cones were biased towards green. These patterns were similar to previous reports from single-cone (Sabesan et al. 2016; Schmidt et al. 2018) and large-field studies (DeValois et al. 1997).

To better appreciate the influence of cone type and number of cones targeted on color reports, data was pooled across subjects and grouped according to the type of cone or pair probed. The mean and standard error for each group is shown in the right most plot below. When an individual or pair of M-cones was targeted the average \(gr\) response was greater than zero indicating a bias towards green. In comparison, the average L-cone(s) elicited biases towards red and yellow. Together these cone type specific differences in color reports were consistent with a predictive relationship between cone type and color report, as previously reported (Sabesan et al. 2016; Schmidt et al. 2018). Two cones with the same photopigment tended to elicit slightly more saturated reports than single cone trials. On the other hand, one L- and one M-cone targeted together tended to produce desaturated reports.

(p <- grid.arrange(ind, avg, widths=2.5:1))

TableGrob (1 x 2) "arrange": 2 grobs
  z     cells    name           grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]

ggsave("figures/UAD_combined.pdf", p)

Mosaic parameters do not predict responses

The plot above demonstrates that color reports varied even between cones with the same photopigment (Schmidt et al. 2018). Some L-cones, for instance, elicited highly saturated red sensations, while a majority produce white or desaturated pink reports. We next asked whether the variability in sensations between cones with the same sensitivity could be explained by low level features of the mosaic. Specifically, can we predict whether an L-cone will produce a saturated red or a desaturated pink based on the surrounding cone types or the delivery error? And in the case of cone pairs, did the distance between the two cones influence color appearance? To answer these question, we grouped responses according to the number of cones targeted (one or two), the session in which trial occurred and the specific cone(s) that was targeted. The mean and count was then computed for each cone(s). Only those cones pairs which had at least four good trials were analyzed.

minTrials <- 4
singleConesSummary <- trial_data %>% 
  select(., -typeIDstr) %>%
  filter(., isPair==FALSE) %>% 
  group_by(., subject, session, masterID1) %>%   
  summarise_all(., c("mean", "length")) %>%
  filter(., yb_length >= minTrials) %>%
  select(., -contains("_length")) %>%
  rename_(.dots=setNames(names(.), gsub("_mean", "", names(.))))
twoConesSummary <- trial_data %>%   
  select(., -typeIDstr) %>%
  filter(., isPair==TRUE) %>% 
  group_by(., subject, session, masterID1, masterID2) %>%
  summarise_all(., c("mean", "length")) %>%
  filter(., yb_length >= minTrials) %>%
  select(., -contains("_length")) %>%
  rename_(.dots=setNames(names(.), gsub("_mean", "", names(.))))

Below is a matrix of correlation plots for the single cone condition.

ptest <- function(mat, ...) {
  # stolen from: http://www.sthda.com/english/wiki/visualize-correlation-matrix-using-correlogram
    mat <- as.matrix(mat)
    n <- ncol(mat)
    p.mat<- matrix(NA, n, n)
    diag(p.mat) <- 0
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            tmp <- cor.test(mat[, i], mat[, j], ...)
            p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
        }
    }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}
corr_analysis <- function (data, predictors) {
  # scatterplotMatrix with car package
  scatterplotMatrix(data[ , predictors], 
                    diagonal=list(method ="histogram"),
                    ellipse=FALSE)
  
  M<-cor(data[ , predictors], use="complete.obs")
  head(round(M,2))
  cor.mtest <- ptest
  # matrix of the p-value of the correlation
  p.mat <- cor.mtest(data[ , predictors])
  print(format(p.mat, digits=4))
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(M, method="color", col=col(200),  
         type="upper", 
         #order="hclust", 
         addCoef.col = "black", # Add coefficient of correlation
         tl.col="black", tl.srt=45, #Text label color and rotation
         # Combine with significance
         p.mat = p.mat, sig.level = 0.01, insig = "blank", 
         # hide correlation coefficient on the principal diagonal
         diag=FALSE 
         )
}
predictors = c("yb", "gr", "saturation", "lConeNeighbors", 
               "distance_to_Scone_1", "distance_to_Scone_2")
corr_analysis(singleConesSummary, predictors)

                    yb          gr          saturation  lConeNeighbors distance_to_Scone_1 distance_to_Scone_2
yb                  "0.000e+00" "2.820e-26" "1.138e-03" "2.000e-01"    "5.135e-01"         "5.398e-01"        
gr                  "2.820e-26" "0.000e+00" "1.101e-01" "2.578e-01"    "9.136e-01"         "9.673e-01"        
saturation          "1.138e-03" "1.101e-01" "0.000e+00" "1.255e-01"    "2.386e-01"         "9.307e-01"        
lConeNeighbors      "2.000e-01" "2.578e-01" "1.255e-01" "0.000e+00"    "1.122e-04"         "1.244e-03"        
distance_to_Scone_1 "5.135e-01" "9.136e-01" "2.386e-01" "1.122e-04"    "0.000e+00"         "4.281e-18"        
distance_to_Scone_2 "5.398e-01" "9.673e-01" "9.307e-01" "1.244e-03"    "4.281e-18"         "0.000e+00"        

Two cone condition:

predictors = c("yb", "gr", "saturation", "lConeNeighbors",
               "distance_btwn_cones_arcmin",  "distance_to_Scone_1", 
               "distance_to_Scone_2")
corr_analysis(twoConesSummary, predictors)

                           yb          gr          saturation  lConeNeighbors distance_btwn_cones_arcmin
yb                         "0.000e+00" "1.963e-31" "9.787e-04" "3.704e-01"    "3.019e-01"               
gr                         "1.963e-31" "0.000e+00" "1.655e-01" "2.613e-02"    "8.605e-01"               
saturation                 "9.787e-04" "1.655e-01" "0.000e+00" "9.544e-01"    "7.302e-03"               
lConeNeighbors             "3.704e-01" "2.613e-02" "9.544e-01" "0.000e+00"    "9.895e-01"               
distance_btwn_cones_arcmin "3.019e-01" "8.605e-01" "7.302e-03" "9.895e-01"    "0.000e+00"               
distance_to_Scone_1        "7.549e-02" "6.574e-01" "7.046e-02" "2.741e-05"    "5.186e-01"               
distance_to_Scone_2        "2.468e-01" "4.834e-01" "6.474e-01" "3.262e-04"    "1.137e-01"               
                           distance_to_Scone_1 distance_to_Scone_2
yb                         "7.549e-02"         "2.468e-01"        
gr                         "6.574e-01"         "4.834e-01"        
saturation                 "7.046e-02"         "6.474e-01"        
lConeNeighbors             "2.741e-05"         "3.262e-04"        
distance_btwn_cones_arcmin "5.186e-01"         "1.137e-01"        
distance_to_Scone_1        "0.000e+00"         "7.036e-17"        
distance_to_Scone_2        "7.036e-17"         "0.000e+00"        

None of the above correlations were statistically significant with the exception of relationships that were known to be dependent, such as between gr and yb. The local neighborhood surrounding a cone is typically thought to be an important factor in generating color sensations. However, we did not find a statistical relationship. The distance between cones may also be an important factor influencing appearance. Neither color nor saturation judgments were correlated with the distance between targeted cones (plotted here in pixels). Cone pairs were never separated by more than one cone, which may explain why we did not detect a relationship. Moreover, the subjects verbally reported that the flashes always appeared as a single uniformly colored dot. Over these very small spatial distances, the visual system appears to be yoking neuronal activities together. In the future, systematically varying the distance between stimulated pairs will be an informative exercise. At a certain critical distance, the points of light will be seen as two spatially distinct dots. It is less clear at what distance the points will be perceived as two distinct colors.

Two cone responses are an average of individual reports

While features of the mosaic and physical stimulus did not predict color reports, we hypothesized that the sensations recorded from individual cones would be predictive of paired stimulation conditions. To address this question, we matched the mean response from each cone pair with the the mean report from each cone tested individually.

merged1 <- merge(twoConesSummary, 
           select(singleConesSummary, subject, 
                  session, masterID1, yb, gr, lConeNeighbors), 
           by=c("subject", "session", "masterID1")) %>% 
  rename(., yb12=yb.x, gr12=gr.x, yb1=yb.y, gr1=gr.y,
                  lConeNeighbors12=lConeNeighbors.x, 
                  lConeNeighbors1=lConeNeighbors.y)
sessionMerge <- merged1 %>% 
  merge(., select(singleConesSummary, subject, session,
                  masterID1, yb, gr, lConeNeighbors), 
        by.x=c("subject", "session", "masterID2"),
        by.y=c("subject", "session", "masterID1")) %>% 
  rename(., yb2=yb, gr2=gr, lConeNeighbors2=lConeNeighbors) %>%
  mutate(., saturation12=abs(yb12) + abs(gr12),
         saturation1=abs(yb1) + abs(gr1),
         saturation2=abs(yb2) + abs(gr2))

We then fit a linear model to the data. Behavioral reports from two cone stimulation were predicted by an average of the individual responses: \(gr_{12} = (gr_1 + gr_2) / 2\). Below is a visualization of the measured responses against the predictions.

sessionMerge['predictedYB'] = (sessionMerge$yb1 + sessionMerge$yb2) / 2
sessionMerge['predictedGR'] = (sessionMerge$gr1 + sessionMerge$gr2) / 2
sessionMerge['predictedSaturation'] = abs(sessionMerge$predictedYB) + abs(sessionMerge$predictedGR)
# Add a column w difference from prediction
sessionMerge <- mutate(sessionMerge, 
                       diff_from_prediction=saturation12 - predictedSaturation)

unity = data.frame(x=c(-1, 1), y=c(-1, 1))
gr_plot <- sessionMerge %>%
  ggplot(., aes(x=predictedGR, y=gr12)) + 
  geom_point() +
  geom_smooth(method=lm) +
  coord_equal(xlim=c(-1, 1), ylim=c(-1, 1)) +
  geom_path(data=unity, aes(x=x, y=y), colour='gray', linetype=1) +
  ggtitle('(green - red) / total') +
  geom_path(data=unity, aes(x=x, y=y), colour='gray', linetype=1) +
  xlab('predicted response') + ylab('') + 
  theme_classic(base_size=15)
unity = data.frame(x=c(-1, 1), y=c(-1, 1))
yb_plot <- sessionMerge %>%
  ggplot(., aes(x=predictedYB, y=yb12)) + 
  geom_point() +
  geom_smooth(method=lm) +
  coord_equal(xlim=c(-1, 1), ylim=c(-1, 1)) +
  geom_path(data=unity, aes(x=x, y=y), colour='gray', linetype=1) +
  ggtitle('(yellow - blue) / total') +
  geom_path(data=unity, aes(x=x, y=y), colour='gray', linetype=1) +
  xlab('predicted response') + 
  ylab('observed pair response') + 
  theme_classic(base_size=15)
p <- grid.arrange(yb_plot, gr_plot, nrow=1)

ggsave("figures/prediction.pdf", p)

Results of linear regression (gr dimension followed by yb):

# A linear model predicts the results based on an average of the individual responses.
mod <- lm(gr12 ~ predictedGR, sessionMerge)
summary.lm(mod)

Call:
lm(formula = gr12 ~ predictedGR, data = sessionMerge)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.45553 -0.09987 -0.00895  0.09816  0.57855 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.01772    0.01259   1.407    0.161    
predictedGR  1.12390    0.05041  22.297   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1623 on 181 degrees of freedom
Multiple R-squared:  0.7331,    Adjusted R-squared:  0.7316 
F-statistic: 497.2 on 1 and 181 DF,  p-value: < 2.2e-16

mod <- lm(yb12 ~ predictedYB, sessionMerge)
summary.lm(mod)

Call:
lm(formula = yb12 ~ predictedYB, data = sessionMerge)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.203925 -0.017762  0.004318  0.018450  0.236417 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.004318   0.004214  -1.025    0.307    
predictedYB  1.002010   0.043201  23.194   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05699 on 181 degrees of freedom
Multiple R-squared:  0.7482,    Adjusted R-squared:  0.7469 
F-statistic:   538 on 1 and 181 DF,  p-value: < 2.2e-16

Pairs with the same type are more saturated than predicted by average

A simple linear model captured a large fraction of the variance (>72%). However, there were some pairs that deviated substantially from the best fit line. We wondered whether the deviation from linearity might be predicted by the sub-class of the two cones. For instance, do an L- and M-cone interact in a non-linear manner, while two L or two M-cones sum linearly?

We found the saturation for each pair and subtracted it from the average of the two cones probed alone. Those results are plotted below. A unity line represents the condition where the observed saturation judgment was predicted exactly by the average of individual responses. Notice that the L+L and M+M pairs tend to fall above the unity line – particularly at higher saturation values. In contrast, the L+M pairs fall below the line. These observations indicate that the cones with the same spectral type produced slightly more saturated reports than predicted by the average of their individual responses.

We quantified this trend directly by taking the difference between the observed and predicted saturation judgments. The results are illustrated in a histogram below. Student’s t-test’s confirm that the L+M pairs were significantly less saturated (more white) than a simple average of their individual responses, while the opposite was true for cones with the same spectral type.

diffHist <- sessionMerge %>%
   filter(., typeID == 1 | typeID > 3) %>%
   ggplot(., 
          aes(x=saturation12 - predictedSaturation, 
              color=as.factor(typeID))) + 
  geom_freqpoly(binwidth=0.1, size=1.25) +
  scale_color_manual(
    values = c(gray, green, yellow, red), 
    name="spectral type",
    labels = c("unknown-pair", "M+M", "L+M", "L+L")) + 
  xlab("measured - predicted saturation") +
  theme_classic(base_size=15) +
  theme(legend.position=c(0.75, 0.80))
ggsave('figures/observed-predicted_histogram.pdf')

TableGrob (1 x 2) "arrange": 2 grobs
  z     cells    name           grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]

L+M pairs had a mean that was significantly different from zero.

t.test(sessionMerge %>% 
         filter(., typeID == 5) %>% 
         select(., diff_from_prediction)
       )

    One Sample t-test

data:  sessionMerge %>% filter(., typeID == 5) %>% select(., diff_from_prediction)
t = -1.9437, df = 56, p-value = 0.05696
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.068670699  0.001035897
sample estimates:
 mean of x 
-0.0338174 

L+L and M+M pairs had a mean that was significantly different from zero, but in the opposite direction.

t.test(sessionMerge %>% 
         filter(., typeID == 4 | typeID == 6) %>%
         select(., diff_from_prediction)
       )

    One Sample t-test

data:  sessionMerge %>% filter(., typeID == 4 | typeID == 6) %>% select(.,     diff_from_prediction)
t = 4.172, df = 78, p-value = 7.772e-05
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.03748812 0.10592157
sample estimates:
 mean of x 
0.07170484 

Unclassified cones were significantly more saturated than predicted.

t.test(sessionMerge %>% 
         filter(., typeID == 1) %>% 
         select(., diff_from_prediction)
       )

    One Sample t-test

data:  sessionMerge %>% filter(., typeID == 1) %>% select(., diff_from_prediction)
t = 4.0435, df = 46, p-value = 0.0001988
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.06116617 0.18243140
sample estimates:
mean of x 
0.1217988 

Finally, L+M pairs were significantly different from the population of responses collected from L+L and M+M pairs.

t.test(sessionMerge %>% 
         filter(., typeID == 5) %>% 
         select(., diff_from_prediction),
       sessionMerge %>% 
         filter(., typeID == 4 | typeID == 6) %>%
         select(., diff_from_prediction)
       )

    Welch Two Sample t-test

data:  sessionMerge %>% filter(., typeID == 5) %>% select(., diff_from_prediction) and sessionMerge %>% filter(., typeID == 4 | typeID == 6) %>% select(., sessionMerge %>% filter(., typeID == 5) %>% select(., diff_from_prediction) and     diff_from_prediction)
t = -4.3148, df = 129.85, p-value = 3.139e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.15390621 -0.05713828
sample estimates:
  mean of x   mean of y 
-0.03381740  0.07170484 

Conclusions

Color sensations from pairs of cones were predicted by a simple average of individual responses. However, when two cones from the same subclass were probed, there was a systematic deviation from a simple average. These pairs produced significantly more saturated colors than predicted by an average of the colors elicited when probed alone. This observation suggests that the visual system uses a different strategy when combining information within versus across neuronal sub-classes.