OSCR

Binocular combination in the autonomic nervous system.

Code ↔ Paper

16 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 16 matches
  1. [1] § Materials and methods › Apparatus and stimuli ↔ manuscript.Rmd, lines 2265–2299 · score 1.00 · depth optical cylinders, circular neutral density, appropriately threaded adapters, guarantee safe illumination, white glass discs, highly diffusive white
  2. [2] § Materials and methods › Apparatus and stimuli ↔ manuscript.Rmd, lines 2265–2299 · score 1.00 · elicit visible percepts, penumbral cone activation, whilst keeping activity, nominally silenced photoreceptors, whilst keeping cone, intrinsically photoreceptive retinal
  3. [3] § Materials and methods › Apparatus and stimuli ↔ manuscript.Rmd, lines 2235–2263 · score 1.00 · Ocean Optics Jaz, NIST traceable light, degree cone fundamentals, distinct spectral tuning, Silent substitution exploits, Argon lamp
  4. [4] § Materials and methods › Computational model and parameter estimation ↔ manuscript.Rmd, lines 2301–2361 · score 1.00 · Target eye assignment, consensual response matched, hierarchical Bayesian framework, Gaussian priors, suppression parameters reached, noise parameter
  5. [5] § Results ↔ manuscript.Rmd, lines 2409–2415 · score 0.99 · downstream temporal filtering, close phase correspondence, rod activity relatively, cone signals differed, quadrature phase, phase arrangement
  6. [6] § Results ↔ manuscript.Rmd, lines 2378–2397 · score 0.99 · Bayesian hierarchical implementation, weakest suppression corresponds, posterior predictive check, fitted model parameters, compare posterior distributions, strongest suppressive weight
  7. [7] § Materials and methods › Apparatus and stimuli ↔ manuscript.Rmd, lines 2235–2263 · score 0.99 · supports synchronous launch, customised binocular viewing, synchronised stimulus modulations, spectral sequences prepared, LEDMOTIVE Technologies, active cooling
  8. [8] § Materials and methods › Procedure ↔ manuscript.Rmd, lines 2301–2361 · score 0.99 · unmodulated background luminance, modified VR headset, black fixation dot, windowless room, 15–17 minutes, 45–60 minutes
  9. [9] § Results ↔ manuscript.Rmd, lines 2378–2397 · score 0.98 · yellow diamond symbols, strongest dichoptic masking, dichoptic cross response, progressive decline, binocular signal combination, strong suppression reducing
  10. [10] § Results ↔ manuscript.Rmd, lines 2399–2407 · score 0.97 · periodic response magnitude, root sum square, extreme ratios, Binocular stimulation exceeded, harmonics moderated, qualitative pattern
  11. [11] § Results ↔ manuscript.Rmd, lines 2364–2376 · score 0.94 · Fourier amplitude spectrum, pupil waveforms, axis scale, light flux modulation, harmonic response, strong modulation
  12. [12] § Results ↔ manuscript.Rmd, lines 2364–2376 · score 0.92 · purely linear system, nonlinear temporal processing, separate oscillatory mechanism, sinusoid responds, nonlinear transformation, harmonic distortion
  13. [13] § Results ↔ manuscript.Rmd, lines 2399–2407 · score 0.81 · weakest suppressive weight, opposite pattern, substantial overlap, lower rows, strongest binocular facilitation, posterior distributions
  14. [14] § Materials and methods › Apparatus and stimuli ↔ manuscript.Rmd, lines 1091–1126 · score 0.61 · opic irradiances, left eye, right eye, power, spectrum, photoreceptor
  15. [15] § Materials and methods › Computational model and parameter estimation ↔ manuscript.Rmd, lines 2417–2452 · score 0.56 · stronger binocular facilitation, ocularity invariance, strong suppression, interocular suppression, binocular response, strength
  16. [16] § Materials and methods › Data analysis ↔ manuscript.Rmd, lines 26–70 · score 0.51 · Mahalanobis distance, analysed, bootstrapping, Fourier, model

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R Markdown · 2,467 lines · 143 KB · no license · 16 matches

  1. ---
  2. title: 'Binocular combination in the autonomic nervous system'
  3. author: |
  4. | Federico G. Segala$^1$*, Aurelio Bruno$^{1,2}$, Joel T. Martin$^{1,3}$,
  5. | Anisa Y. Morsi$^1$, Alex R. Wade$^{1,4}$ & Daniel H. Baker$^{1,4}$
  6. output:
  7. bookdown::pdf_document2:
  8. fig_caption: yes
  9. toc: no
  10. keep_tex: true
  11. html_document: default
  12. bibliography: references.bib
  13. link-citations: yes
  14. linkcolor: black
  15. csl: plos-one.csl
  16. header-includes:
  17. \usepackage{caption}
  18. \usepackage{setspace}\doublespacing
  19. \usepackage{lineno} \linenumbers
  20. \newcommand{\beginsupplement}{\setcounter{table}{0}
  21. \renewcommand{\thetable}{S\arabic{table}}
  22. \setcounter{figure}{0}
  23. \renewcommand{\thefigure}{S\arabic{figure}}}
  24. ---
  26. ```{r setup, include=FALSE}
  28. processdata <- 0 # this flag determines the amount of processing, with 4 levels:
  29. # 0 - do no processing, generate the pdf using existing versions of all figures
  30. # 1 - generate figures using the processed group data
  31. # 2 - average data across participants, do bootstrapping
  32. # 3 - download all raw data and analyse (requires 56GB of storage)
  34. domodelling <- 0 # if set to 1, run the Stan models (takes ~24 hours)
  36. # all files are stored in the project directory /local/ which git is told to ignore
  37. # create a local directory to store data and outputs
  38. localdir <- 'local/'
  39. if (!file.exists(localdir)){dir.create(localdir)}
  41. # create a local directory to store raw data
  42. rawdir <- 'local/rawdata/'
  43. if (!file.exists(rawdir)){dir.create(rawdir)}
  45. # create a local directory to store processed data
  46. datadir <- 'local/processeddata/'
  47. if (!file.exists(datadir)){dir.create(datadir)}
  49. # create a local directory to store figures
  50. figdir <- 'Figures/'
  51. if (!file.exists(figdir)){dir.create(figdir)}
  53. # cut off for data inclusion (Mahalanobis distance for complex data)
  54. SDthresh <<- 3
  55. nbootstraps <<- 10000
  57. # check which packages are installed, install the missing ones, and activate
  58. packagelist <- c('knitr','remotes','tictoc','R.matlab','bookdown','grImport','png','tiff','pals','ez','gtools','signal','boot','rstan','bayestestR','rstatix','coda','parallel','utils','osfr','grid','gridExtra','cowplot','ggsignif') # list of CRAN packages
  60. missingpackages <- packagelist[!packagelist %in% installed.packages()[,1]]
  61. if (length(missingpackages)>0){install.packages(missingpackages,repos = "http://cran.us.r-project.org")}
  62. if (!'FourierStats' %in% installed.packages()[,1]){remotes::install_github("bakerdh/FourierStats")}
  63. packagelist <- c(packagelist,'FourierStats')
  64. toinstall <- packagelist[which(!packagelist %in% (.packages()))]
  65. invisible(lapply(toinstall,library,character.only=TRUE))
  67. devtools::install_github("kupietz/kableExtra")
  68. library(kableExtra)
  70. knitr::opts_chunk$set(echo = TRUE)
  72. # settings for sampler
  73. numSavedSteps <- 100000 # total number of MCMC steps
  74. adaptSteps <- 1000 # Number of steps to "tune" the samplers
  75. burnInSteps <- 2000
  76. thinSteps <- 10
  77. nChains <- 16
  78. ncores <- 8
  79. options(mc.cores=ncores)
  81. # helper function to make colours transparent
  82. addalpha <- function(col, alpha=1){apply(sapply(col, col2rgb)/255, 2, function(x) rgb(x[1], x[2], x[3], alpha=alpha))}
  84. anovacirc2way <- function(data, groupA=NULL, groupB=NULL, participant=NULL){
  86. grouplabelsA <- groupA
  87. grouplabelsB <- groupB
  88. participantlabels <- participant
  89. datavals <- data
  91. grouplabelsA <- as.factor(grouplabelsA)
  92. factorlistA <- levels(grouplabelsA)
  93. grouplabelsB <- as.factor(grouplabelsB)
  94. factorlistB <- levels(grouplabelsB)
  95. participantlabels <- as.factor(participantlabels)
  96. participantlist <- levels(participantlabels)
  98. grandmean <- mean(datavals)
  99. grandvar <- sum(abs(datavals - grandmean)^2)/(length(datavals)-1)
  101. SST <- grandvar*(length(datavals) - 1)
  103. SSW <- 0
  104. for (n in 1:nlevels(participantlabels)){SSW <- SSW + sum(abs(datavals[which(participantlabels==participantlist[n])] - mean(datavals[which(participantlabels==participantlist[n])]))^2)}
  105. dfW <- 2*(nlevels(participantlabels)*(nlevels(grouplabelsA)*nlevels(grouplabelsB)-1))
  107. SSBet <- SST - SSW # between subjects SS - not actually used
  108. SSM <- 0
  109. for (a in 1:nlevels(grouplabelsA)){
  110. for (b in 1:nlevels(grouplabelsB)){
  111. i <- as.numeric(((grouplabelsA==grouplabelsA[a])))
  112. j <- as.numeric(((grouplabelsB==grouplabelsB[b])))
  113. indices <- which((i*j)>0)
  114. SSM <- SSM + (nlevels(participantlabels)*sum(abs(datavals[indices] - grandmean)^2))
  115. }
  116. }
  117. dfM <- 2*((nlevels(grouplabelsA)*nlevels(grouplabelsB))-1)
  119. SSR <- 0
  120. for (a in 1:nlevels(grouplabelsA)){
  121. for (b in 1:nlevels(grouplabelsB)){
  122. i <- as.numeric(((grouplabelsA==grouplabelsA[a])))
  123. j <- as.numeric(((grouplabelsB==grouplabelsB[b])))
  124. indices <- which((i*j)>0)
  125. condvar <- sum(abs(datavals[indices] - mean(datavals[indices]))^2)/(length(indices)-1)
  126. SSR <- SSR + (nlevels(participantlabels)-1)*condvar
  127. }
  128. }
  129. dfR <- 2*(nlevels(participantlabels)-1)*(nlevels(grouplabelsA)*nlevels(grouplabelsB))
  131. groupmeansA <- NULL
  132. for (n in 1:nlevels(grouplabelsA)){groupmeansA[n] <- mean(datavals[which(grouplabelsA==factorlistA[n])])}
  133. SSA <- 0
  134. for (n in 1:nlevels(grouplabelsA)){SSA <- SSA + length(which(grouplabelsA==factorlistA[n]))*abs(groupmeansA[n]-grandmean)^2}
  135. dfA <- 2*(nlevels(grouplabelsA)-1)
  137. groupmeansB <- NULL
  138. for (n in 1:nlevels(grouplabelsB)){groupmeansB[n] <- mean(datavals[which(grouplabelsB==factorlistB[n])])}
  139. SSB <- 0
  140. for (n in 1:nlevels(grouplabelsB)){SSB <- SSB + length(which(grouplabelsB==factorlistB[n]))*abs(groupmeansB[n]-grandmean)^2}
  141. dfB <- 2*(nlevels(grouplabelsB)-1)
  143. SSAB <- SSW - (SSA + SSB)
  144. dfAB <- 2 * (nlevels(grouplabelsA)-1) * (nlevels(grouplabelsB)-1)
  146. MSA <- SSA/dfA
  147. MSB <- SSB/dfB
  148. MSAB <- SSAB/dfAB
  149. MSR <- SSR/dfR
  151. FratioA <- MSA/MSR
  152. pvalA <- pf(FratioA,df1=dfA,df2=dfR,lower.tail=FALSE)
  153. FratioB <- MSB/MSR
  154. pvalB <- pf(FratioB,df1=dfB,df2=dfR,lower.tail=FALSE)
  155. FratioAB <- MSAB/MSR
  156. pvalAB <- pf(FratioAB,df1=dfAB,df2=dfR,lower.tail=FALSE)
  158. method <- 'Factorial repeated measures ANOVA^2-circ'
  159. output <- data.frame(FratioA,pvalA,FratioB,pvalB,FratioAB,pvalAB,dfA,dfB,dfAB,dfR,method)
  161. return(output)
  162. }
  164. cleanup2D <- function(data){
  165. # helper function to clean up complex data by removing any data points more than
  166. # SDthresh standard deviations from the mean (in Mahalanobis distance units)
  167. cartdata <- data.frame(Re(data),Im(data))
  168. output <- 0
  169. if (sum(cartdata)!=0){
  170. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  171. i <- which(D<SDthresh)
  172. output <- mean(temp[i])}
  173. return(output)}
  175. bootstrap2D <- function(data){
  176. mbs <- data.frame(meanamp=0,lowerCI=0,upperCI=0)
  177. cartdata <- data.frame(Re(data),Im(data))
  178. if (sum(cartdata)!=0){
  179. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  180. i <- which(D<SDthresh)
  181. mbs <- amperrors(data[i],method='boot',quantiles=68,nresamples=nbootstraps)}
  182. return(mbs)}
  184. contrastsdB <- 20*log10(c(6,12,24,48,96))
  185. colpal <- c('#FE5000','#8783CF','#228B22','#808080','#6d008b','#8B8000')
  186. colpalmod <- c('#808080','#8B8000','#FE5000','#8783CF')
  188. # these are approximate values for isolumiant red, green, blue and yellow
  189. stimcolours <- c(rgb(0.02,0.02,0.02),rgb(0.98,0.98,0.98),rgb(0,0.61,0.5),rgb(1,0.39,0.5),rgb(0.35,0.63,0),rgb(0.65,0.37,1),'grey','cornflowerblue')
  191. if (processdata>2){
  192. # check if raw data are present locally, and download if not
  193. nodelist <- c('6u7db','xkgd9','pexyk','k9mdq') # list of OSF nodes containing raw data
  194. for (expt in 1:4){
  195. osfproject <- osf_retrieve_node(nodelist[expt])
  196. osffiles <- osf_ls_files(osfproject,n_max=300)
  197. for (f in 1:nrow(osffiles)){
  198. fname <- osffiles$name[f]
  199. foldername <- strsplit(fname,'.',fixed=TRUE)[[1]][1]
  200. if (!dir.exists(paste0(rawdir,foldername))){
  201. osf_download(osffiles[f,],localdir,progress=TRUE)
  202. dir.create(paste0(rawdir,foldername))
  203. unzip(paste0(localdir,fname),exdir=paste0(rawdir,foldername))
  204. file.remove(paste0(localdir,fname))
  205. }
  206. }
  207. }
  208. }
  210. if (processdata==2 || domodelling==1){
  211. # if we're going to need the pre-processed files, extract from zip file
  212. if (!file.exists(paste0(datadir,'P301_summary.RData'))){
  213. if (!file.exists(paste0(localdir,'IndividualSummary.zip'))){
  214. osfnode <- '9ebtd'
  215. osfproject <- osf_retrieve_node(osfnode)
  216. osffiles <- osf_ls_files(osfproject,n_max=300)
  217. fid <- which(osffiles$name=='IndividualSummary.zip')
  218. osf_download(osffiles[fid,],localdir,progress=TRUE)
  219. }
  220. unzip(paste0(localdir,'IndividualSummary.zip'),exdir=datadir)
  221. }
  222. }
  224. if (processdata<2){
  225. # if we're using the pre-processed data, check it's present and download if not
  226. exptnames <- c('AveragedataLum.RData','AveragedataLM.RData','AveragedataS.RData','AveragedataMel.RData','included.RData')
  227. for (f in 1:length(exptnames)){
  228. if (!file.exists(paste0(datadir,exptnames[f]))){
  229. if (!exists('osffiles')){
  230. osfnode <- '9ebtd'
  231. osfproject <- osf_retrieve_node(osfnode)
  232. osffiles <- osf_ls_files(osfproject,n_max=300)
  233. }
  234. fid <- which(osffiles$name==exptnames[f])
  235. osf_download(osffiles[fid,],datadir,progress=TRUE)
  236. }}}
  238. if (domodelling==0){
  239. # if we're not doing the modelling, check if the model outputs are present and download if not
  240. exptnames <- c('modeloutputLum.RData','modeloutputLM.RData','modeloutputS.RData','modeloutputMel.RData')
  242. for (f in 1:length(exptnames)){
  243. if (!file.exists(paste0(datadir,exptnames[f]))){
  244. if (!exists('osffiles')){
  245. osfnode <- '9ebtd'
  246. osfproject <- osf_retrieve_node(osfnode)
  247. osffiles <- osf_ls_files(osfproject,n_max=300)
  248. }
  249. fid <- which(osffiles$name==exptnames[f])
  250. osf_download(osffiles[fid,],datadir,progress=TRUE)
  251. }
  252. }
  253. }
  255. ```
  257. $^1$Department of Psychology, University of York, Heslington, York, YO10 5DD, United Kingdom.
  259. $^2$School of Psychology and Vision Sciences, University of Leicester, Leicester, LE1 7RH, United Kingdom.
  261. $^3$School of Philosophy, Psychology and Language Sciences, University of Edinburgh, Edinburgh, EH8 9AD, United Kingdom.
  263. $^4$York Biomedical Research Institute, University of York, Heslington, York, YO10 5NG, United Kingdom.
  265. *Corresponding author
  267. E-mail:[[email hidden]]([email hidden]) (FGS)
  269. \pagebreak
  271. ```{r individualanalyses, include=FALSE, results='hide'}
  273. # analyse individual participant data for all conditions
  275. if (processdata > 2){
  277. allsubjects <- dir(rawdir)
  278. for (p in 1:length(allsubjects)){
  279. participant <- allsubjects[p]
  281. pupiltargets <- array(0,dim=c(3,2,60))
  282. pupilmasks <- pupiltargets
  283. pupiltargets2 <- pupiltargets
  284. pupilmasks2 <- pupiltargets
  285. pupilwaveforms <- array(0,dim=c(3,2,60,120*14))
  286. pupilspectra <- array(0,dim=c(3,2,60,300))
  287. timeseq <- seq(1/120,10,length.out=120*10)
  288. timeseq2 <- seq(1/120,14,length.out=120*14)
  289. targetindex <- (0.5*10)+1
  290. maskindex <- (0.4*10)+1
  291. targetindex2 <- (0.5*2*10)+1
  292. maskindex2 <- (2*0.4*10)+1
  294. repcounterP <- (1:60)*0
  295. annfiles <- dir(paste0(rawdir,participant),pattern='*annotations.csv',full.names=TRUE)
  296. pupilfiles <- dir(paste0(rawdir,participant),pattern='*pupil_positions.csv',full.names=TRUE)
  298. for (block in 1:length(annfiles)){
  300. adata <- read.csv(annfiles[block])
  302. trialtimes <- NULL
  303. trialtypes <- NULL
  304. for (t in 1:(nrow(adata)/2)){
  305. trialtimes[t] <- mean(adata$timestamp[(t*2 - 1):(t*2)])
  306. trialtypes[t] <- adata$condition_code[t*2]
  307. }
  309. trialtypes[which(trialtypes>15)] <- trialtypes[which(trialtypes>15)] + 5
  310. inclist <- NULL
  311. for (t in 11:15){
  312. i <- which(trialtypes==t)
  313. inclist[t-10] <- i[2]
  314. }
  315. inclist <- inclist[which(!is.na(inclist))]
  316. trialtypes[inclist] <- trialtypes[inclist] + 5
  318. pdata <- read.csv(pupilfiles[block])
  319. pdata2 <- pdata[,c(1,3,4,14)]
  320. pdata2 <- pdata2[which(pdata2[,3]>0.2),]
  322. for (trialno in 1:length(trialtypes)){
  323. repcounterP[trialtypes[trialno]] <- repcounterP[trialtypes[trialno]] + 1
  324. for (eye in 1:2){
  325. eyedata <- pdata2[which(pdata2[,2]==(2-eye)),]
  326. a <- which(eyedata[,1]>trialtimes[trialno]+2)
  327. b <- which(eyedata[,1]<(trialtimes[trialno]+12))
  328. i <- intersect(a,b)
  329. trial <- eyedata[i,]
  330. trial <- trial[which(!is.na(trial[,4])),]
  331. trial[,1] <- trial[,1] - trial[1,1]
  332. if (nrow(trial)>3){
  333. resampled <- interp1(trial[,1],trial[,4],timeseq,method='linear',extrap=TRUE)
  334. fspec <- (fft(resampled)/length(resampled))
  335. pupiltargets[repcounterP[trialtypes[trialno]],eye,trialtypes[trialno]] <- fspec[targetindex]
  336. pupilmasks[repcounterP[trialtypes[trialno]],eye,trialtypes[trialno]] <- fspec[maskindex]
  337. pupiltargets2[repcounterP[trialtypes[trialno]],eye,trialtypes[trialno]] <- fspec[targetindex2]
  338. pupilmasks2[repcounterP[trialtypes[trialno]],eye,trialtypes[trialno]] <- fspec[maskindex2]
  339. pupilspectra[repcounterP[trialtypes[trialno]],eye,trialtypes[trialno],1:300] <- fspec[1:300]
  340. }
  342. a <- which(eyedata[,1]>(trialtimes[trialno]-1))
  343. b <- which(eyedata[,1]<(trialtimes[trialno]+13))
  344. i <- intersect(a,b)
  345. trial <- eyedata[i,]
  346. trial[,1] <- trial[,1] - trial[1,1]
  347. if (nrow(trial)>3){
  348. trial <- trial[which(!is.na(trial[,4])),]
  349. resampled <- interp1(trial[,1],trial[,4],timeseq2,method='linear',extrap=TRUE)
  350. pupilwaveforms[repcounterP[trialtypes[trialno]],eye,trialtypes[trialno],] <- resampled - mean(resampled[1:120])
  351. }
  352. }
  353. }
  354. }
  356. cleanmeansP <- matrix(0,nrow=6,ncol=5)
  357. cleanmasksP <- matrix(0,nrow=6,ncol=5)
  358. cleanmeansP2 <- matrix(0,nrow=6,ncol=5)
  359. cleanmasksP2 <- matrix(0,nrow=6,ncol=5)
  361. for (cond in 1:6){
  362. startindex <- (10*(cond-1))
  363. for (level in 1:5){
  364. temp <- c(pupiltargets[,,startindex+level],pupiltargets[,,startindex+level+5])
  365. cartdata <- data.frame(Re(temp),Im(temp))
  366. if (sum(cartdata)!=0){
  367. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  368. i <- which(D<SDthresh)
  369. cleanmeansP[cond,level] <- mean(temp[i])}
  371. temp <- c(pupiltargets2[,,startindex+level],pupiltargets2[,,startindex+level+5])
  372. cartdata <- data.frame(Re(temp),Im(temp))
  373. if (sum(cartdata)!=0){
  374. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  375. i <- which(D<SDthresh)
  376. cleanmeansP2[cond,level] <- mean(temp[i])}
  378. temp <- c(pupilmasks[,,startindex+level],pupilmasks[,,startindex+level+5])
  379. cartdata <- data.frame(Re(temp),Im(temp))
  380. if (sum(cartdata)!=0){
  381. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  382. i <- which(D<SDthresh)
  383. cleanmasksP[cond,level] <- mean(temp[i])}
  385. temp <- c(pupilmasks2[,,startindex+level],pupilmasks2[,,startindex+level+5])
  386. cartdata <- data.frame(Re(temp),Im(temp))
  387. if (sum(cartdata)!=0){
  388. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  389. i <- which(D<SDthresh)
  390. cleanmasksP2[cond,level] <- mean(temp[i])}
  391. }
  392. }
  394. meanwavesP <- apply(pupilwaveforms[,,c(15,20),],4,mean,na.rm=TRUE)
  395. meanspectraP <- apply(pupilspectra[,,c(15,20),],4,mean,na.rm=TRUE)
  396. save(file=paste(datadir,participant,'_summary.RData',sep=''),list=c('cleanmeansP','cleanmasksP','cleanmeansP2','cleanmasksP2','meanspectraP','meanwavesP'))
  397. }
  399. }
  401. ```
  403. ```{r groupanalyses, include=FALSE, results='hide'}
  405. # files from P301-P324 are for S-cone stimulation
  406. # files from P501-P524 are for melanopsin stimulation
  407. # files from P601-P624 are for luminance stimulation
  408. # files from P701-P724 are for (L-M) cone stimulation
  410. # load in individual participant data and average, save group data
  411. if (processdata > 1){
  413. prefixlist <- c('P6','P7','P3','P5')
  414. exptnames <- c('AveragedataLum.RData','AveragedataLM.RData','AveragedataS.RData','AveragedataMel.RData')
  416. includedvalues <- array(0,dim=c(4,24,6,5))
  417. includedvalues2 <- array(0,dim=c(4,24,6,5))
  418. includedmasks <- array(0,dim=c(4,24,6,5))
  419. includedmasks2 <- array(0,dim=c(4,24,6,5))
  421. for (expt in 1:4){
  422. d <- dir(datadir,pattern=prefixlist[expt], full.names = TRUE)
  424. tempmasksP <- array(0,dim=c(length(d),6,5))
  425. tempmeansP <- tempmasksP
  426. tempmeansP2 <- tempmasksP
  427. tempmasksP2 <- tempmasksP
  428. tempspectraP <- array(0,dim=c(length(d),300))
  429. tempwavesP <- array(0,dim=c(length(d),1680))
  431. for (s in 1:length(d)){
  433. load(d[s])
  435. tempmasksP[s,,] <- cleanmasksP
  436. tempmeansP[s,,] <- cleanmeansP
  437. tempmasksP2[s,,] <- cleanmasksP2
  438. tempmeansP2[s,,] <- cleanmeansP2
  439. tempspectraP[s,] <- meanspectraP
  440. tempwavesP[s,] <- meanwavesP
  441. }
  443. # tempspectraP <- tempspectraP[-c(2,4,13,20,26),]
  445. cleanmeansP <- matrix(0,nrow=6,ncol=5)
  446. cleanmeansPCI <- array(0,c(2,6,5))
  447. cleanmasksP <- matrix(0,nrow=6,ncol=5)
  448. cleanmasksPCI <- array(0,c(2,6,5))
  449. cleanmeansP2 <- matrix(0,nrow=6,ncol=5)
  450. cleanmeansP2CI <- array(0,c(2,6,5))
  451. cleanmasksP2 <- matrix(0,nrow=6,ncol=5)
  452. cleanmasksP2CI <- array(0,c(2,6,5))
  454. for (cond in 1:6){
  455. for (level in 1:5){
  457. temp <- tempmeansP[,cond,level]
  458. cartdata <- data.frame(Re(temp),Im(temp))
  459. if (sum(cartdata)!=0){
  460. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  461. i <- which(D<SDthresh)
  462. includedvalues[expt,i,cond,level] <- 1
  463. mbs <- amperrors(temp[i],method='boot',quantiles=68,nresamples=nbootstraps)
  464. cleanmeansP[cond,level] <- mean(temp[i])
  465. cleanmeansPCI[1,cond,level] <- mbs$lowerCI
  466. cleanmeansPCI[2,cond,level] <- mbs$upperCI
  467. }
  469. temp <- tempmasksP[,cond,level]
  470. cartdata <- data.frame(Re(temp),Im(temp))
  471. if (sum(cartdata)!=0){
  472. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  473. i <- which(D<SDthresh)
  474. includedmasks[expt,i,cond,level] <- 1
  475. mbs <- amperrors(temp[i],method='boot',quantiles=68,nresamples=nbootstraps)
  476. cleanmasksP[cond,level] <- mean(temp[i])
  477. cleanmasksPCI[1,cond,level] <- mbs$lowerCI
  478. cleanmasksPCI[2,cond,level] <- mbs$upperCI
  479. }
  481. temp <- tempmeansP2[,cond,level]
  482. cartdata <- data.frame(Re(temp),Im(temp))
  483. if (sum(cartdata)!=0){
  484. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  485. i <- which(D<SDthresh)
  486. includedvalues2[expt,i,cond,level] <- 1
  487. mbs <- amperrors(temp[i],method='boot',quantiles=68,nresamples=nbootstraps)
  488. cleanmeansP2[cond,level] <- mean(temp[i])
  489. cleanmeansP2CI[1,cond,level] <- mbs$lowerCI
  490. cleanmeansP2CI[2,cond,level] <- mbs$upperCI
  491. }
  493. temp <- tempmasksP2[,cond,level]
  494. cartdata <- data.frame(Re(temp),Im(temp))
  495. if (sum(cartdata)!=0){
  496. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  497. i <- which(D<SDthresh)
  498. includedmasks2[expt,i,cond,level] <- 1
  499. mbs <- amperrors(temp[i],method='boot',quantiles=68,nresamples=nbootstraps)
  500. cleanmasksP2[cond,level] <- mean(temp[i])
  501. cleanmasksP2CI[1,cond,level] <- mbs$lowerCI
  502. cleanmasksP2CI[2,cond,level] <- mbs$upperCI
  503. }
  504. }
  505. }
  507. meanspectraP <- (1:300)*0
  508. meanspectraPCI <- array(0,c(2,300))
  510. for (f in 2:300){
  511. temp <- tempspectraP[,f]
  512. cartdata <- data.frame(Re(temp),Im(temp))
  513. if (sum(cartdata)!=0){
  514. D <- sqrt(stats::mahalanobis(cartdata, colMeans(cartdata), cov(cartdata)))
  515. i <- which(D<SDthresh)
  516. mbs <- amperrors(temp[i],method='boot',quantiles=68,nresamples=nbootstraps)
  517. meanspectraP[f] <- mbs$meanamp
  518. meanspectraPCI[1,f] <- mbs$lowerCI
  519. meanspectraPCI[2,f] <- mbs$upperCI
  520. }
  521. }
  523. meanwavesP <- NULL
  524. meanwavesPCI <- array(0,dim=c(2,dim(tempwavesP)[2]))
  526. nsubjs <- dim(tempwavesP)[1]
  527. for (t in 1:dim(tempwavesP)[2]){
  528. temp <- tempwavesP[,t]
  529. i <- which(abs(temp-mean(temp))<(SDthresh*sd(temp)))
  530. meanwavesP[t] <- mean(temp[i])
  531. bspop <- NULL
  532. for (s in 1:nbootstraps){bspop[s] <- mean(sample(temp[i],length(i),replace=TRUE))}
  533. meanwavesPCI[,t] <- quantile(bspop,c(0.16,0.84),na.rm=TRUE)
  534. }
  536. save(file=paste0(datadir,exptnames[expt]),list=c('meanwavesP','meanwavesPCI','cleanmeansP','cleanmeansPCI','cleanmasksP','cleanmasksPCI','cleanmeansP2','cleanmeansP2CI','cleanmasksP2','cleanmasksP2CI','meanspectraP','meanspectraPCI','tempmeansP','tempmeansP2','tempmasksP','tempmasksP2'))
  537. }
  539. save(file=paste0(datadir,'included.RData'),list=c('includedvalues','includedvalues2','includedmasks','includedmasks2'))
  541. }
  542. ```
  544. ```{r makestanmodels, include=FALSE, results='hide'}
  546. # create text objects for Stan models
  548. if (domodelling > 0){
  550. CRFmodelstring = "
  551. data {
  552. int<lower=1> Nsubj ;
  553. int<lower=1> Ntotal ;
  554. real y[Ntotal] ;
  555. real x[Ntotal] ;
  556. int<lower=1> c[Ntotal] ;
  557. int<lower=1> s[Ntotal] ;
  558. }
  559. parameters {
  560. real<lower=1> Z[Nsubj] ;
  561. real<lower=0> Rmax[Nsubj] ;
  562. real<lower=0> wt[Nsubj] ;
  563. real<lower=0> k[Nsubj] ;
  564. real<lower=0> p[Nsubj] ;
  565. real<lower=0> q[Nsubj] ;
  566. real<lower=0> sigma ;
  567. real<lower=1> Zmu ;
  568. real<lower=0> Rmaxmu ;
  569. real<lower=0> wtmu ;
  570. real<lower=0> kmu ;
  571. real<lower=0> pmu ;
  572. real<lower=0> qmu ;
  573. real<lower=0> Zsigma ;
  574. real<lower=0> Rmaxsigma ;
  575. real<lower=0> wtsigma ;
  576. real<lower=0> ksigma ;
  577. real<lower=0> psigma ;
  578. real<lower=0> qsigma ;
  579. real<lower=0> nu ;
  580. }
  581. model {
  582. Zmu ~ normal( 40 , 20 ) ;
  583. Rmaxmu ~ normal( 0.5 , 0.3 ) ;
  584. wtmu ~ normal( 1 , 0.5 ) ;
  585. kmu ~ normal( 0.002 , 0.001 ) ;
  586. pmu ~ normal( 0.4 , 0.2 ) ;
  587. qmu ~ normal( 2 , 0.5 ) ;
  588. sigma ~ uniform( 1.0E-3 , 1.0E+3 ) ;
  589. Zsigma ~ uniform( 1.0E-3 , 1.0E+3 ) ;
  590. Rmaxsigma ~ uniform( 1.0E-3 , 1.0E+3 ) ;
  591. wtsigma ~ uniform( 1.0E-3 , 1.0E+3 ) ;
  592. ksigma ~ uniform( 1.0E-3 , 1.0E+3 ) ;
  593. psigma ~ uniform( 1.0E-3 , 1.0E+3 ) ;
  594. qsigma ~ uniform( 1.0E-3 , 1.0E+3 ) ;
  595. nu ~ exponential(1/30.0) ;
  596. Z ~ normal( Zmu , Zsigma ) ; // vectorized
  597. Rmax ~ normal( Rmaxmu , Rmaxsigma ) ; // vectorized
  598. wt ~ normal( wtmu , wtsigma ) ; // vectorized
  599. k ~ normal( kmu , ksigma ) ; // vectorized
  600. p ~ normal( pmu, psigma ) ;
  601. q ~ normal( qmu, qsigma ) ;
  602. for ( i in 1:Ntotal ) {
  604. if (c[i]==1) // monocular condition F1
  605. y[i] ~ student_t(
  606. nu ,
  607. k[s[i]] + (Rmax[s[i]]*(pow(x[i],p[s[i]]+q[s[i]]))/(Z[s[i]] + (pow(x[i],q[s[i]])))) ,
  608. sigma ) ;
  610. else if (c[i]==2) // binocular condition F1
  611. y[i] ~ student_t(
  612. nu ,
  613. k[s[i]] + (Rmax[s[i]]*((pow(x[i],p[s[i]]+q[s[i]]))/(Z[s[i]] + (pow(x[i],q[s[i]])) + wt[s[i]]*(pow(x[i],q[s[i]]))) + (pow(x[i],p[s[i]]+q[s[i]]))/(Z[s[i]] + (pow(x[i],q[s[i]])) + wt[s[i]]*(pow(x[i],q[s[i]]))))),
  614. sigma ) ;
  616. else if (c[i]==3) // dichoptic condition F1
  617. y[i] ~ student_t(
  618. nu ,
  619. k[s[i]] + (Rmax[s[i]]*((pow(x[i],p[s[i]]+q[s[i]]))/(Z[s[i]] + (pow(x[i],q[s[i]])) + wt[s[i]]*(pow(48,q[s[i]]))) + (pow(48,p[s[i]]+q[s[i]]))/(Z[s[i]] + (pow(48,q[s[i]])) + wt[s[i]]*(pow(x[i],q[s[i]]))))),
  620. sigma ) ;
  623. else if (c[i]==4) // binocular cross condition F1
  624. y[i] ~ student_t(
  625. nu ,
  626. k[s[i]] + (Rmax[s[i]]*(pow(x[i],p[s[i]]+q[s[i]]))/(Z[s[i]] + (pow(x[i],q[s[i]])) + wt[s[i]]*(pow(x[i],q[s[i]])))) ,
  627. sigma ) ;
  630. else if (c[i]==5) // dichoptic cross condition F1
  631. y[i] ~ student_t(
  632. nu ,
  633. k[s[i]] + (Rmax[s[i]]*(pow(x[i],p[s[i]]+q[s[i]]))/(Z[s[i]] + (pow(x[i],q[s[i]])) + wt[s[i]]*(pow(48,q[s[i]])))) ,
  634. sigma ) ;
  637. else if (c[i]==6) // monocular condition F2
  638. y[i] ~ student_t(
  639. nu ,
  640. k[s[i]] + (Rmax[s[i]]*(pow(x[i],p[s[i]]+q[s[i]]))/(Z[s[i]] + (pow(x[i],q[s[i]])))) ,
  641. sigma ) ;
  643. else if (c[i]==7) // binocular cross condition F2
  644. y[i] ~ student_t(
  645. nu ,
  646. k[s[i]] + (Rmax[s[i]]*(pow(x[i],p[s[i]]+q[s[i]]))/(Z[s[i]] + (pow(x[i],q[s[i]])) + wt[s[i]]*(pow(x[i],q[s[i]])))) ,
  647. sigma ) ;
  650. else if (c[i]==8) // dichoptic cross condition F2
  651. y[i] ~ student_t(
  652. nu ,
  653. k[s[i]] + (Rmax[s[i]]*(pow(48,p[s[i]]+q[s[i]]))/(Z[s[i]] + (pow(48,q[s[i]])) + wt[s[i]]*(pow(x[i],q[s[i]])))) ,
  654. sigma ) ;
  656. }
  657. }
  659. " # close quote for CRFmodelstring
  661. # Translate Stan models to C++ and compile to DSO:
  662. CRFmodel <- stan_model(model_code=CRFmodelstring)
  663. }
  666. ```
  668. ```{r runstanmodels, include=FALSE, results='hide'}
  670. # do computational modelling for all experiments
  672. # files from P301-P324 are for S-cone stimulation
  673. # files from P501-P524 are for melanopsin stimulation
  674. # files from P601-P624 are for luminance stimulation
  675. # files from P701-P724 are for (L-M) cone stimulation
  677. if (domodelling > 0){
  679. prefixlist <- c('P6','P7','P3','P5')
  680. exptnames <- c('modeloutputLum.RData','modeloutputLM.RData','modeloutputS.RData','modeloutputMel.RData')
  681. load(paste0(datadir,'included.RData'))
  683. for (expt in 1:length(prefixlist)){
  685. d <- dir(datadir,pattern=prefixlist[expt])
  687. allmasksP <- array(0,dim=c(length(d),6,5))
  688. allmeansP <- allmasksP
  689. allmasksP2 <- array(0,dim=c(length(d),6,5))
  690. allmeansP2 <- allmasksP
  692. for (s in 1:length(d)){
  694. load(paste0(datadir,d[s]))
  696. allmasksP[s,,] <- cleanmasksP
  697. allmeansP[s,,] <- cleanmeansP
  698. allmasksP2[s,,] <- cleanmasksP2
  699. allmeansP2[s,,] <- cleanmeansP2
  701. }
  703. # 1st harmonic F1
  704. condlist <- c(1,2,3,5,6)
  705. contlist <- c(6,12,24,48,96)
  706. y <- NULL
  707. x <- NULL
  708. c <- NULL
  709. s <- NULL
  710. i <- 0
  711. for (subj in 1:dim(allmeansP)[1]){
  712. for (cond in 1:length(condlist)){
  713. for (cont in 1:length(contlist)){
  714. dp <- abs(allmeansP[subj,condlist[cond],cont])
  715. if (includedvalues[expt,subj,condlist[cond],cont]==1){
  716. if (is.na(dp)==0){
  717. if (is.infinite(dp)==0){
  718. i <- i + 1
  719. x[i] <- contlist[cont]
  720. y[i] <- dp
  721. c[i] <- cond
  722. s[i] <- subj
  723. }
  724. }
  725. }
  726. }
  727. }
  728. }
  729. # 1st harmonic F2
  730. condlist <- c(4,5,6)
  731. for (subj in 1:dim(allmasksP)[1]){
  732. for (cond in 1:length(condlist)){
  733. for (cont in 1:length(contlist)){
  734. dp <- abs(allmasksP[subj,condlist[cond],cont])
  735. if (includedmasks[expt,subj,condlist[cond],cont]==1){
  736. if (is.na(dp)==0){
  737. if (is.infinite(dp)==0){
  738. i <- i + 1
  739. x[i] <- contlist[cont]
  740. y[i] <- dp
  741. c[i] <- cond+5
  742. s[i] <- subj
  743. }
  744. }
  745. }
  746. }
  747. }
  748. }
  752. pupildatalist = list(x = x, y = y, s = s, c = c, Nsubj = max(s), Ntotal = length(y))
  754. # Get MC sample of posterior:
  755. pupilsamples <- sampling(object=CRFmodel,
  756. data = pupildatalist,
  757. chains = nChains,
  758. iter = (ceiling(numSavedSteps/nChains)*thinSteps+burnInSteps),
  759. warmup = burnInSteps,
  760. thin = thinSteps,
  761. cores = getOption("mc.cores", 1L))
  763. # 2nd harmonic
  764. condlist <- c(1,2,3,5,6)
  765. contlist <- c(6,12,24,48,96)
  766. y <- NULL
  767. x <- NULL
  768. c <- NULL
  769. s <- NULL
  770. i <- 0
  771. for (subj in 1:dim(allmeansP2)[1]){
  772. for (cond in 1:length(condlist)){
  773. for (cont in 1:length(contlist)){
  774. dp <- abs(allmeansP2[subj,condlist[cond],cont])
  775. if (includedvalues2[expt,subj,condlist[cond],cont]==1){
  776. if (is.na(dp)==0){
  777. if (is.infinite(dp)==0){
  778. i <- i + 1
  779. x[i] <- contlist[cont]
  780. y[i] <- dp
  781. c[i] <- cond
  782. s[i] <- subj
  783. }
  784. }
  785. }
  786. }
  787. }
  788. }
  789. # 2nd harmonic F2
  790. condlist <- c(4,5,6)
  791. for (subj in 1:dim(allmasksP)[1]){
  792. for (cond in 1:length(condlist)){
  793. for (cont in 1:length(contlist)){
  794. dp <- abs(allmasksP2[subj,condlist[cond],cont])
  795. if (includedmasks2[expt,subj,condlist[cond],cont]==1){
  796. if (is.na(dp)==0){
  797. if (is.infinite(dp)==0){
  798. i <- i + 1
  799. x[i] <- contlist[cont]
  800. y[i] <- dp
  801. c[i] <- cond+5
  802. s[i] <- subj
  803. }
  804. }
  805. }
  806. }
  807. }
  808. }
  810. pupildatalist = list(x = x, y = y, s = s, c = c, Nsubj = max(s), Ntotal = length(y))
  812. # Get MC sample of posterior:
  813. pupilsamples2 <- sampling(object=CRFmodel,
  814. data = pupildatalist,
  815. chains = nChains,
  816. iter = (ceiling(numSavedSteps/nChains)*thinSteps+burnInSteps),
  817. warmup = burnInSteps,
  818. thin = thinSteps,
  819. cores = getOption("mc.cores", 1L))
  821. save(file=paste0(datadir,exptnames[expt]),list=c('pupilsamples','pupilsamples2'))
  822. }
  824. }
  826. ```
  828. ```{r getmodelpredictions, include=FALSE, results='hide'}
  830. # get predictions of the model based on earlier fits
  832. load(paste0(datadir,'modeloutputLum.RData'))
  834. draws <- extract(pupilsamples)
  836. pupilweightpopLum <- draws$wtmu
  837. pupilZpopLum <- draws$Zmu
  838. pupilKpopLum <- draws$kmu
  839. pupilPpopLum <- draws$pmu
  840. pupilQpopLum <- draws$qmu
  841. pupilRpopLum <- draws$Rmaxmu
  843. pupilWLum <- 10^as.numeric(map_estimate(log10(pupilweightpopLum)))
  844. pupilZLum <- 10^as.numeric(map_estimate(log10(pupilZpopLum)))
  845. pupilKLum <- 10^as.numeric(map_estimate(log10(pupilKpopLum)))/sqrt(24)
  846. pupilPLum <- 10^as.numeric(map_estimate(log10(pupilPpopLum+pupilQpopLum)))
  847. pupilQLum <- 10^as.numeric(map_estimate(log10(pupilQpopLum)))
  848. pupilRmaxLum <- 10^as.numeric(map_estimate(log10(pupilRpopLum)))
  850. load(paste0(datadir,'modeloutputMel.RData'))
  852. draws <- extract(pupilsamples)
  854. pupilweightpopMel <- draws$wtmu
  855. pupilZpopMel <- draws$Zmu
  856. pupilKpopMel <- draws$kmu
  857. pupilPpopMel <- draws$pmu
  858. pupilQpopMel <- draws$qmu
  859. pupilRpopMel <- draws$Rmaxmu
  861. pupilWMel <- 10^as.numeric(map_estimate(log10(pupilweightpopMel)))
  862. pupilZMel <- 10^as.numeric(map_estimate(log10(pupilZpopMel)))
  863. pupilKMel <- 10^as.numeric(map_estimate(log10(pupilKpopMel)))/sqrt(24)
  864. pupilPMel <- 10^as.numeric(map_estimate(log10(pupilPpopMel+pupilQpopMel)))
  865. pupilQMel <- 10^as.numeric(map_estimate(log10(pupilQpopMel)))
  866. pupilRmaxMel <- 10^as.numeric(map_estimate(log10(pupilRpopMel)))
  868. load(paste0(datadir,'modeloutputLM.RData'))
  870. draws <- extract(pupilsamples)
  872. pupilweightpopLM <- draws$wtmu
  873. pupilZpopLM <- draws$Zmu
  874. pupilKpopLM <- draws$kmu
  875. pupilPpopLM <- draws$pmu
  876. pupilQpopLM <- draws$qmu
  877. pupilRpopLM <- draws$Rmaxmu
  879. pupilWLM <- 10^as.numeric(map_estimate(log10(pupilweightpopLM)))
  880. pupilZLM <- 10^as.numeric(map_estimate(log10(pupilZpopLM)))
  881. pupilKLM <- 10^as.numeric(map_estimate(log10(pupilKpopLM)))/sqrt(24)
  882. pupilPLM <- 10^as.numeric(map_estimate(log10(pupilPpopLM+pupilQpopLM)))
  883. pupilQLM <- 10^as.numeric(map_estimate(log10(pupilQpopLM)))
  884. pupilRmaxLM <- 10^as.numeric(map_estimate(log10(pupilRpopLM)))
  886. load(paste0(datadir,'modeloutputS.RData'))
  888. draws <- extract(pupilsamples)
  890. pupilweightpopS <- draws$wtmu
  891. pupilZpopS <- draws$Zmu
  892. pupilKpopS <- draws$kmu
  893. pupilPpopS <- draws$pmu
  894. pupilQpopS <- draws$qmu
  895. pupilRpopS <- draws$Rmaxmu
  897. pupilWS <- 10^as.numeric(map_estimate(log10(pupilweightpopS)))
  898. pupilZS <- 10^as.numeric(map_estimate(log10(pupilZpopS)))
  899. pupilKS <- 10^as.numeric(map_estimate(log10(pupilKpopS)))/sqrt(24)
  900. pupilPS <- 10^as.numeric(map_estimate(log10(pupilPpopS+pupilQpopS)))
  901. pupilQS <- 10^as.numeric(map_estimate(log10(pupilQpopS)))
  902. pupilRmaxS <- 10^as.numeric(map_estimate(log10(pupilRpopS)))
  904. load(paste0(datadir,'modeloutputLum.RData'))
  906. draws <- extract(pupilsamples2)
  908. pupilweightpopLum2 <- draws$wtmu
  909. pupilZpopLum2 <- draws$Zmu
  910. pupilKpopLum2 <- draws$kmu
  911. pupilPpopLum2 <- draws$pmu
  912. pupilQpopLum2 <- draws$qmu
  913. pupilRpopLum2 <- draws$Rmaxmu
  915. pupilWLum2 <- 10^as.numeric(map_estimate(log10(pupilweightpopLum2)))
  916. pupilZLum2 <- 10^as.numeric(map_estimate(log10(pupilZpopLum2)))
  917. pupilKLum2 <- 10^as.numeric(map_estimate(log10(pupilKpopLum2)))/sqrt(24)
  918. pupilPLum2 <- 10^as.numeric(map_estimate(log10(pupilPpopLum2+pupilQpopLum2)))
  919. pupilQLum2 <- 10^as.numeric(map_estimate(log10(pupilQpopLum2)))
  920. pupilRmaxLum2 <- 10^as.numeric(map_estimate(log10(pupilRpopLum2)))
  922. load(paste0(datadir,'modeloutputMel.RData'))
  924. draws <- extract(pupilsamples2)
  926. pupilweightpopMel2 <- draws$wtmu
  927. pupilZpopMel2 <- draws$Zmu
  928. pupilKpopMel2 <- draws$kmu
  929. pupilPpopMel2 <- draws$pmu
  930. pupilQpopMel2 <- draws$qmu
  931. pupilRpopMel2 <- draws$Rmaxmu
  933. pupilWMel2 <- 10^as.numeric(map_estimate(log10(pupilweightpopMel2)))
  934. pupilZMel2 <- 10^as.numeric(map_estimate(log10(pupilZpopMel2)))
  935. pupilKMel2 <- 10^as.numeric(map_estimate(log10(pupilKpopMel2)))/sqrt(24)
  936. pupilPMel2 <- 10^as.numeric(map_estimate(log10(pupilPpopMel2+pupilQpopMel2)))
  937. pupilQMel2 <- 10^as.numeric(map_estimate(log10(pupilQpopMel2)))
  938. pupilRmaxMel2 <- 10^as.numeric(map_estimate(log10(pupilRpopMel2)))
  940. load(paste0(datadir,'modeloutputLM.RData'))
  942. draws <- extract(pupilsamples2)
  944. pupilweightpopLM2 <- draws$wtmu
  945. pupilZpopLM2 <- draws$Zmu
  946. pupilKpopLM2 <- draws$kmu
  947. pupilPpopLM2 <- draws$pmu
  948. pupilQpopLM2 <- draws$qmu
  949. pupilRpopLM2 <- draws$Rmaxmu
  951. pupilWLM2 <- 10^as.numeric(map_estimate(log10(pupilweightpopLM2)))
  952. pupilZLM2 <- 10^as.numeric(map_estimate(log10(pupilZpopLM2)))
  953. pupilKLM2 <- 10^as.numeric(map_estimate(log10(pupilKpopLM2)))/sqrt(24)
  954. pupilPLM2 <- 10^as.numeric(map_estimate(log10(pupilPpopLM2+pupilQpopLM2)))
  955. pupilQLM2 <- 10^as.numeric(map_estimate(log10(pupilQpopLM2)))
  956. pupilRmaxLM2 <- 10^as.numeric(map_estimate(log10(pupilRpopLM2)))
  958. load(paste0(datadir,'modeloutputS.RData'))
  960. draws <- extract(pupilsamples2)
  962. pupilweightpopS2 <- draws$wtmu
  963. pupilZpopS2 <- draws$Zmu
  964. pupilKpopS2 <- draws$kmu
  965. pupilPpopS2 <- draws$pmu
  966. pupilQpopS2 <- draws$qmu
  967. pupilRpopS2 <- draws$Rmaxmu
  969. pupilWS2 <- 10^as.numeric(map_estimate(log10(pupilweightpopS2)))
  970. pupilZS2 <- 10^as.numeric(map_estimate(log10(pupilZpopS2)))
  971. pupilKS2 <- 10^as.numeric(map_estimate(log10(pupilKpopS2)))/sqrt(24)
  972. pupilPS2 <- 10^as.numeric(map_estimate(log10(pupilPpopS2+pupilQpopS2)))
  973. pupilQS2 <- 10^as.numeric(map_estimate(log10(pupilQpopS2)))
  974. pupilRmaxS2 <- 10^as.numeric(map_estimate(log10(pupilRpopS2)))
  976. getmodresps <- function(params){
  977. contrastsfinedB <- 2:40
  978. contrastsfine <- 10^(contrastsfinedB/20)
  980. modresps <- matrix(0,nrow=5,ncol=length(contrastsfine))
  981. modresps[1,] <- params[1]*(contrastsfine^params[2])/(params[3] + contrastsfine^params[6]) + params[5]
  983. binresp <- 2*((contrastsfine^params[2])/(params[3] + contrastsfine^params[6] + params[4]*contrastsfine^params[6]))
  984. modresps[2,] <- params[1]*binresp + params[5]
  986. dichresp <- ((contrastsfine^params[2])/(params[3] + contrastsfine^params[6] + params[4]*48^params[6])) + ((48^params[2])/(params[3] + 48^params[6] + params[4]*contrastsfine^params[6]))
  987. modresps[3,] <- params[1]*dichresp + params[5]
  989. dichresp2 <- ((contrastsfine^params[2])/(params[3] + contrastsfine^params[6] + params[4]*48^params[6]))
  990. modresps[4,] <- params[1]*dichresp2 + params[5]
  992. dichresp3 <- ((48^params[2])/(params[3] + 48^params[6] + params[4]*contrastsfine^params[6]))
  993. modresps[5,] <- params[1]*dichresp3 + params[5]
  995. return(modresps)}
  997. params <- c(pupilRmaxLum,pupilPLum,pupilZLum,pupilWLum,pupilKLum,pupilQLum)
  998. modrespsLum <- getmodresps(params)
  1000. params <- c(pupilRmaxMel,pupilPMel,pupilZMel,pupilWMel,pupilKMel,pupilQMel)
  1001. modrespsMel <- getmodresps(params)
  1003. params <- c(pupilRmaxLM,pupilPLM,pupilZLM,pupilWLM,pupilKLM,pupilQLM)
  1004. modrespsLM <- getmodresps(params)
  1006. params <- c(pupilRmaxS,pupilPS,pupilZS,pupilWS,pupilKS,pupilQS)
  1007. modrespsS <- getmodresps(params)
  1009. params <- c(pupilRmaxLum2,pupilPLum2,pupilZLum2,pupilWLum2,pupilKLum2,pupilQLum2)
  1010. modrespsLum2 <- getmodresps(params)
  1012. params <- c(pupilRmaxMel2,pupilPMel2,pupilZMel2,pupilWMel2,pupilKMel2,pupilQMel2)
  1013. modrespsMel2 <- getmodresps(params)
  1015. params <- c(pupilRmaxLM2,pupilPLM2,pupilZLM2,pupilWLM2,pupilKLM2,pupilQLM2)
  1016. modrespsLM2 <- getmodresps(params)
  1018. params <- c(pupilRmaxS2,pupilPS2,pupilZS2,pupilWS2,pupilKS2,pupilQS2)
  1019. modrespsS2 <- getmodresps(params)
  1021. ```
  1023. ```{r makemethodsfig, include=FALSE, results='hide'}
  1025. # create plots of spectra, cone absorption, waveforms and Fourier spectra
  1026. if (processdata > 0){
  1028. plotcontrasts <- 1 # either plot photoreceptor contrasts (1) or absolute irradiance change (0)
  1030. times <- seq((1/120)-1,13,1/120)
  1031. frequencies <- (0:299)/10
  1032. sinewave <- sin(0.5*times * 2*pi)
  1033. sinewave[1:120] <- 0
  1034. sinewave[1560:1680] <- 0
  1036. coneab <- read.csv('resources/coneabsorptions.csv',header=FALSE)
  1037. spectra <- read.csv('resources/spectra.csv',header=FALSE)
  1039. # calculate delta photoreceptor contrast
  1040. pcontrast <- coneab
  1041. for (n in 4:11){pcontrast[,n] <- 100*coneab[,n]/coneab[,2]}
  1043. pdf(paste0(figdir,"methodsfigrods.pdf"), bg="transparent", height = 13, width = 12)
  1045. layout(matrix(1:20, 5, 4, byrow = TRUE))
  1046. par(mar=c(1.5,1,0.1,0.1))
  1048. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=c(0,1), ylim=c(0,1))
  1050. schematic <- readTIFF('resources/schematic.tiff')
  1051. aspratio <- 15/12 # this is the aspect ratio of the output pdf
  1052. imwidth <- 0.85
  1053. xstart <- 0
  1054. ystart <- 0
  1055. rasterImage(schematic,xstart,ystart,xstart+imwidth*aspratio,ystart+imwidth)
  1057. text(0,0.88,'(a)',pos=4,cex=2)
  1059. par(mar=c(3.5,6,1.5,2))
  1061. lspec <- read.csv('resources/STLAB_1_oo_irrad_spectra.csv',header=FALSE)
  1062. rspec <- read.csv('resources/STLAB_2_oo_irrad_spectra.csv',header=FALSE)
  1064. lall <- matrix(0,nrow=11,ncol=401)
  1065. lall[1,] <- unlist(lspec[1,3:403])
  1066. rall <- matrix(0,nrow=11,ncol=401)
  1067. rall[1,] <- unlist(rspec[1,3:403])
  1068. for (primary in 1:10){
  1069. temp <- subset(lspec,V1==(primary-1) & V2=='4095')
  1070. lall[primary+1,] <- unlist(temp[1,3:403])
  1071. temp <- subset(rspec,V1==(primary-1) & V2=='4095')
  1072. rall[primary+1,] <- unlist(temp[1,3:403])
  1073. }
  1075. primaryrgbvals <- c(rgb(128/255,50/255,218/255),rgb(74/255,104/255,218/255),rgb(0,0,133/255),rgb(0,0,245/255),rgb(117/255,252/255,253/255),rgb(55/255,126/255,33/255),rgb(117/255,252/255,76/255),rgb(242/255,169/255,59/255),rgb(234/255,51/255,35/255),rgb(127/255,23/255,14/255))
  1077. plotlims <- c(380,800,0,0.2)
  1078. ticklocsx <- seq(400,800,100) # locations of tick marks on x axis
  1079. ticklocsy <- seq(0,0.2,0.1) # locations of tick marks on y axis
  1080. ticklabelsx <- ticklocsx # set labels for x ticks
  1081. ticklabelsy <- ticklocsy # set labels for y ticks
  1083. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1084. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1085. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1086. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1087. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1088. title(xlab="Wavelength (nm)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1089. title(ylab=expression(paste("Power (",mu,"W/cm"^"2","/nm)")), col.lab=rgb(0,0,0), line=2, cex.lab=1.5)
  1091. # for (primary in 1:10){polygon(c(lall[1,],780),c(lall[primary+1,],0),col=addalpha(primaryrgbvals[primary],0.5),border=NA)}
  1092. for (primary in 1:10){polygon(c(rall[1,],780),c(rall[primary+1,],0),col=addalpha(primaryrgbvals[primary],0.5),border=NA)}
  1094. legend(680,0.2,1:10,col=primaryrgbvals,lwd=2,box.lwd=2,title='Primary')
  1095. text(380,0.19,'(b)',pos=4,cex=2)
  1098. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1099. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1100. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1101. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1102. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1103. title(xlab="Wavelength (nm)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1104. title(ylab=expression(paste("Power (",mu,"W/cm"^"2","/nm)")), col.lab=rgb(0,0,0), line=2, cex.lab=1.5)
  1106. lines(spectra[,1],spectra[,2],lwd=2,col=colpal[2])
  1107. lines(spectra[,1],spectra[,3],lwd=2,lty=2)
  1109. legend(380,0.2,c('Left eye', 'Right eye'),lty=1:2,lwd=2,col=c(colpal[2],'black'),box.lwd=2)
  1111. text(720,0.19,'(c)',pos=4,cex=2)
  1114. plotlims <- c(0.5,5.5,0,15)
  1115. ticklocsx <- 1:5 # locations of tick marks on x axis
  1116. ticklocsy <- seq(0,15,5) # locations of tick marks on y axis
  1117. ticklabelsx <- c('sc','mc','lc','rods','mel') # set labels for x ticks
  1118. ticklabelsy <- ticklocsy # set labels for y ticks
  1120. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1121. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1122. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1123. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1124. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1125. title(xlab="Photoreceptor class", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1126. title(ylab=expression(paste(alpha,"-opic irradiance")), col.lab=rgb(0,0,0), line=1.8, cex.lab=1.5)
  1128. # coneab <- coneab[c(1:3,5),] # don't plot rod predictions
  1129. for (b in 1:5){
  1130. polygon(c(b-0.3,b,b,b-0.3),c(0,0,coneab[b,2],coneab[b,2]),col=colpal[2],border=NA)
  1131. polygon(c(b+0.3,b,b,b+0.3),c(0,0,coneab[b,3],coneab[b,3]),col='black',border=NA)
  1132. }
  1134. text(4.5,14,'(d)',pos=4,cex=2)
  1137. plotlims <- c(380,800,-0.2,0.2)
  1138. ticklocsx <- seq(400,800,100) # locations of tick marks on x axis
  1139. ticklocsy <- seq(-0.2,0.2,0.1) # locations of tick marks on y axis
  1140. ticklabelsx <- ticklocsx # set labels for x ticks
  1141. ticklabelsy <- ticklocsy # set labels for y ticks
  1143. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1144. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1145. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1146. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1147. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1148. title(xlab="Wavelength (nm)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1149. title(ylab="Difference power", col.lab=rgb(0,0,0), line=2.8, cex.lab=1.5)
  1150. title(main='Light flux',cex.main=2)
  1152. polygon(c(spectra[,1],spectra[401,1]),c(spectra[,4],0),border=NA,col=stimcolours[1])
  1153. polygon(c(spectra[,1],spectra[401,1]),c(spectra[,5],0),border=NA,col=stimcolours[2])
  1155. lines(spectra[,1],spectra[,4],lwd=2)
  1156. lines(spectra[,1],spectra[,5],lwd=2)
  1158. lines(c(380,800),c(0,0),lwd=2)
  1160. text(720,0.18,'(e)',pos=4,cex=2)
  1162. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1163. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1164. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1165. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1166. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1167. title(xlab="Wavelength (nm)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1168. title(ylab="Difference power", col.lab=rgb(0,0,0), line=2.8, cex.lab=1.5)
  1169. title(main='L-M cone',cex.main=2)
  1171. polygon(spectra[,1],spectra[,6],border=NA,col=stimcolours[3])
  1172. polygon(spectra[,1],spectra[,7],border=NA,col=stimcolours[4])
  1174. lines(spectra[,1],spectra[,6],lwd=2)
  1175. lines(spectra[,1],spectra[,7],lwd=2)
  1177. lines(c(380,800),c(0,0),lwd=2)
  1179. text(720,0.18,'(f)',pos=4,cex=2)
  1181. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1182. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1183. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1184. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1185. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1186. title(xlab="Wavelength (nm)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1187. title(ylab="Difference power", col.lab=rgb(0,0,0), line=2.8, cex.lab=1.5)
  1188. title(main='S-cone',cex.main=2)
  1190. polygon(spectra[,1],spectra[,8],border=NA,col=stimcolours[5])
  1191. polygon(spectra[,1],spectra[,9],border=NA,col=stimcolours[6])
  1193. lines(spectra[,1],spectra[,8],lwd=2)
  1194. lines(spectra[,1],spectra[,9],lwd=2)
  1196. lines(c(380,800),c(0,0),lwd=2)
  1198. text(720,0.18,'(g)',pos=4,cex=2)
  1200. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1201. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1202. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1203. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1204. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1205. title(xlab="Wavelength (nm)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1206. title(ylab="Difference power", col.lab=rgb(0,0,0), line=2.8, cex.lab=1.5)
  1207. title(main='Melanopsin',cex.main=2)
  1209. polygon(spectra[,1],spectra[,10],border=NA,col=stimcolours[7])
  1210. polygon(spectra[,1],spectra[,11],border=NA,col=stimcolours[8])
  1212. lines(spectra[,1],spectra[,10],lwd=2)
  1213. lines(spectra[,1],spectra[,11],lwd=2)
  1215. lines(c(380,800),c(0,0),lwd=2)
  1217. text(720,0.18,'(h)',pos=4,cex=2)
  1220. if (plotcontrasts==1){
  1222. # pcontrast <- pcontrast[c(1,2,3,5),] # don't plot rods
  1224. plotlims <- c(0.5,5.5,-100,100)
  1225. ticklocsx <- 1:5 # locations of tick marks on x axis
  1226. ticklocsy <- seq(-100,100,50) # locations of tick marks on y axis
  1227. ticklabelsx <- c('sc','mc','lc','rods','mel') # set labels for x ticks
  1228. ticklabelsy <- ticklocsy # set labels for y ticks
  1230. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1231. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1232. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1233. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1234. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1235. title(xlab="Photoreceptor class", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1236. title(ylab="Contrast (%)", col.lab=rgb(0,0,0), line=2, cex.lab=1.5)
  1238. for (b in 1:5){
  1239. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,pcontrast[b,4],pcontrast[b,4]),col=stimcolours[1],border='black',lwd=2)
  1240. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,pcontrast[b,5],pcontrast[b,5]),col=stimcolours[2],border='black',lwd=2)
  1241. }
  1243. lines(c(0.5,5.5),c(0,0),lwd=2)
  1245. text(4.65,80,'(i)',pos=4,cex=2)
  1247. plotlims <- c(0.5,5.5,-20,20)
  1248. ticklocsx <- 1:5 # locations of tick marks on x axis
  1249. ticklocsy <- seq(-20,20,10) # locations of tick marks on y axis
  1250. ticklabelsy <- ticklocsy # set labels for y ticks
  1252. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1253. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1254. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1255. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1256. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1257. title(xlab="Photoreceptor class", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1258. title(ylab="Contrast (%)", col.lab=rgb(0,0,0), line=2, cex.lab=1.5)
  1260. for (b in 1:5){
  1261. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,pcontrast[b,6],pcontrast[b,6]),col=stimcolours[3],border='black',lwd=2)
  1262. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,pcontrast[b,7],pcontrast[b,7]),col=stimcolours[4],border='black',lwd=2)
  1263. }
  1265. lines(c(0.5,5.5),c(0,0),lwd=2)
  1267. text(4.6,18,'(j)',pos=4,cex=2)
  1269. plotlims <- c(0.5,5.5,-50,50)
  1270. ticklocsx <- 1:5 # locations of tick marks on x axis
  1271. ticklocsy <- seq(-50,50,25) # locations of tick marks on y axis
  1272. ticklabelsy <- ticklocsy # set labels for y ticks
  1274. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1275. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1276. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1277. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1278. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1279. title(xlab="Photoreceptor class", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1280. title(ylab="Contrast (%)", col.lab=rgb(0,0,0), line=2, cex.lab=1.5)
  1282. for (b in 1:5){
  1283. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,pcontrast[b,8],pcontrast[b,8]),col=stimcolours[5],border='black',lwd=2)
  1284. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,pcontrast[b,9],pcontrast[b,9]),col=stimcolours[6],border='black',lwd=2)
  1285. }
  1287. lines(c(0.5,5.5),c(0,0),lwd=2)
  1289. text(4.6,45,'(k)',pos=4,cex=2)
  1291. plotlims <- c(0.5,5.5,-20,20)
  1292. ticklocsx <- 1:5 # locations of tick marks on x axis
  1293. ticklocsy <- seq(-20,20,10) # locations of tick marks on y axis
  1294. ticklabelsy <- ticklocsy # set labels for y ticks
  1296. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1297. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1298. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1299. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1300. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1301. title(xlab="Photoreceptor class", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1302. title(ylab="Contrast (%)", col.lab=rgb(0,0,0), line=2, cex.lab=1.5)
  1304. for (b in 1:5){
  1305. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,pcontrast[b,10],pcontrast[b,10]),col=stimcolours[7],border='black',lwd=2)
  1306. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,pcontrast[b,11],pcontrast[b,11]),col=stimcolours[8],border='black',lwd=2)
  1307. }
  1309. lines(c(0.5,5.5),c(0,0),lwd=2)
  1311. text(0.5,18,'(l)',pos=4,cex=2)
  1312. }
  1314. if (plotcontrasts==0){
  1315. plotlims <- c(0.5,4.5,-15,15)
  1316. ticklocsx <- 1:4 # locations of tick marks on x axis
  1317. ticklocsy <- seq(-15,15,5) # locations of tick marks on y axis
  1318. ticklabelsx <- c('sc','mc','lc','mel') # set labels for x ticks
  1319. ticklabelsy <- ticklocsy # set labels for y ticks
  1321. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1322. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1323. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1324. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1325. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1326. title(xlab="Photoreceptor class", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1327. title(ylab="Irradiance change", col.lab=rgb(0,0,0), line=2.8, cex.lab=1.5)
  1329. for (b in 1:4){
  1330. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,coneab[b,4],coneab[b,4]),col=stimcolours[1],border='black',lwd=2)
  1331. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,coneab[b,5],coneab[b,5]),col=stimcolours[2],border='black',lwd=2)
  1332. }
  1334. lines(c(0.5,4.5),c(0,0),lwd=2)
  1336. text(3.5,14,'(i)',pos=4,cex=2)
  1338. plotlims <- c(0.5,4.5,-4,4)
  1339. ticklocsx <- 1:4 # locations of tick marks on x axis
  1340. ticklocsy <- seq(-4,4,2) # locations of tick marks on y axis
  1341. ticklabelsy <- ticklocsy # set labels for y ticks
  1343. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1344. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1345. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1346. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1347. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1348. title(xlab="Photoreceptor class", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1349. title(ylab="Irradiance change", col.lab=rgb(0,0,0), line=1.8, cex.lab=1.5)
  1351. for (b in 1:4){
  1352. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,coneab[b,6],coneab[b,6]),col=stimcolours[3],border='black',lwd=2)
  1353. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,coneab[b,7],coneab[b,7]),col=stimcolours[4],border='black',lwd=2)
  1354. }
  1356. lines(c(0.5,4.5),c(0,0),lwd=2)
  1358. text(3.5,3.6,'(j)',pos=4,cex=2)
  1360. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1361. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1362. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1363. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1364. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1365. title(xlab="Photoreceptor class", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1366. title(ylab="Irradiance change", col.lab=rgb(0,0,0), line=1.8, cex.lab=1.5)
  1368. for (b in 1:4){
  1369. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,coneab[b,8],coneab[b,8]),col=stimcolours[5],border='black',lwd=2)
  1370. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,coneab[b,9],coneab[b,9]),col=stimcolours[6],border='black',lwd=2)
  1371. }
  1373. lines(c(0.5,4.5),c(0,0),lwd=2)
  1375. text(3.5,3.6,'(k)',pos=4,cex=2)
  1377. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1378. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1379. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1380. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1381. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1382. title(xlab="Photoreceptor class", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1383. title(ylab="Irradiance change", col.lab=rgb(0,0,0), line=1.8, cex.lab=1.5)
  1385. for (b in 1:4){
  1386. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,coneab[b,10],coneab[b,10]),col=stimcolours[7],border='black',lwd=2)
  1387. polygon(c(b-0.3,b+0.3,b+0.3,b-0.3),c(0,0,coneab[b,11],coneab[b,11]),col=stimcolours[8],border='black',lwd=2)
  1388. }
  1390. lines(c(0.5,4.5),c(0,0),lwd=2)
  1392. text(3.5,3.6,'(l)',pos=4,cex=2)
  1393. }
  1396. prefixlist <- c('P6','P7','P3','P5')
  1398. d <- dir(datadir,pattern=prefixlist[1], full.names = TRUE)
  1401. load(file=paste0(datadir,'AveragedataLum.RData'))
  1403. plotlims <- c(-1,13,-1,0.5)
  1404. ticklocsx <- c(-1,seq(0,12,2),13) # locations of tick marks on x axis
  1405. ticklocsy <- seq(-1,0.5,0.5) # locations of tick marks on y axis
  1406. ticklabelsx <- c('',seq(0,12,2),'') # set labels for x ticks
  1407. ticklabelsy <- ticklocsy # set labels for y ticks
  1409. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1410. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1411. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1412. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1413. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1414. title(xlab="Time (s)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1415. title(ylab="Diameter change (mm)", col.lab=rgb(0,0,0), line=2.8, cex.lab=1.5)
  1417. filt <- butter(3,0.083)
  1418. filtwaveP <- signal::filter(filt, meanwavesP)
  1419. filtU <- signal::filter(filt,meanwavesPCI[1,])
  1420. filtL <- signal::filter(filt,meanwavesPCI[2,])
  1421. polygon(times[c(1:length(times),length(times):1)],c(filtL,filtU[length(filtU):1]),col=addalpha(colpal[2],0.4),border=NA)
  1422. lines(times, filtwaveP, col=colpal[2], lwd=1.5)
  1423. lines(times, (sinewave/10)-0.9, col='black', lwd=2)
  1425. for (part in 1:24) {
  1426. load(file=paste0(d[part]))
  1427. filt <- butter(3,0.083)
  1428. filtwaveP <- signal::filter(filt, meanwavesP)
  1429. lines(times, filtwaveP, col=addalpha(colpal[2],0.7), lwd=0.2)
  1430. }
  1432. text(-1,0.4,'(m)',pos=4,cex=2)
  1436. load(file=paste0(datadir,'AveragedataLM.RData'))
  1437. d <- dir(datadir,pattern=prefixlist[2], full.names = TRUE)
  1439. plotlims <- c(-1,13,-1,0.5)
  1440. ticklocsx <- c(-1,seq(0,12,2),13) # locations of tick marks on x axis
  1441. ticklocsy <- seq(-1,0.5,0.5) # locations of tick marks on y axis
  1442. ticklabelsx <- c('',seq(0,12,2),'') # set labels for x ticks
  1443. ticklabelsy <- ticklocsy # set labels for y ticks
  1445. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1446. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1447. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1448. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1449. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1450. title(xlab="Time (s)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1451. title(ylab="Diameter change (mm)", col.lab=rgb(0,0,0), line=2.8, cex.lab=1.5)
  1453. filt <- butter(3,0.083)
  1454. filtwaveP <- signal::filter(filt, meanwavesP)
  1455. filtU <- signal::filter(filt,meanwavesPCI[1,])
  1456. filtL <- signal::filter(filt,meanwavesPCI[2,])
  1457. polygon(times[c(1:length(times),length(times):1)],c(filtL,filtU[length(filtU):1]),col=addalpha(colpal[2],0.4),border=NA)
  1458. lines(times, filtwaveP, col=colpal[2], lwd=1.5)
  1459. lines(times, (sinewave/10)-0.9, col='black', lwd=2)
  1461. for (part in 1:24) {
  1462. load(file=paste0(d[part]))
  1463. filt <- butter(3,0.083)
  1464. filtwaveP <- signal::filter(filt, meanwavesP)
  1465. lines(times, filtwaveP, col=addalpha(colpal[2],0.7), lwd=0.2)
  1466. }
  1468. text(-1,0.4,'(n)',pos=4,cex=2)
  1470. load(file=paste0(datadir,'AveragedataS.RData'))
  1471. d <- dir(datadir,pattern=prefixlist[3], full.names = TRUE)
  1473. plotlims <- c(-1,13,-1,0.5)
  1474. ticklocsx <- c(-1,seq(0,12,2),13) # locations of tick marks on x axis
  1475. ticklocsy <- seq(-1,0.5,0.5) # locations of tick marks on y axis
  1476. ticklabelsx <- c('',seq(0,12,2),'') # set labels for x ticks
  1477. ticklabelsy <- ticklocsy # set labels for y ticks
  1479. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1480. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1481. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1482. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1483. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1484. title(xlab="Time (s)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1485. title(ylab="Diameter change (mm)", col.lab=rgb(0,0,0), line=2.8, cex.lab=1.5)
  1487. filt <- butter(3,0.083)
  1488. filtwaveP <- signal::filter(filt, meanwavesP)
  1489. filtU <- signal::filter(filt,meanwavesPCI[1,])
  1490. filtL <- signal::filter(filt,meanwavesPCI[2,])
  1491. polygon(times[c(1:length(times),length(times):1)],c(filtL,filtU[length(filtU):1]),col=addalpha(colpal[2],0.4),border=NA)
  1492. lines(times, filtwaveP, col=colpal[2], lwd=1.5)
  1493. lines(times, (sinewave/10)-0.9, col='black', lwd=2)
  1495. for (part in 1:24) {
  1496. load(file=paste0(d[part]))
  1497. filt <- butter(3,0.083)
  1498. filtwaveP <- signal::filter(filt, meanwavesP)
  1499. lines(times, filtwaveP, col=addalpha(colpal[2],0.7), lwd=0.2)
  1500. }
  1502. text(-1,0.4,'(o)',pos=4,cex=2)
  1504. load(file=paste0(datadir,'AveragedataMel.RData'))
  1505. d <- dir(datadir,pattern=prefixlist[4], full.names = TRUE)
  1507. plotlims <- c(-1,13,-1,0.5)
  1508. ticklocsx <- c(-1,seq(0,12,2),13) # locations of tick marks on x axis
  1509. ticklocsy <- seq(-1,0.5,0.5) # locations of tick marks on y axis
  1510. ticklabelsx <- c('',seq(0,12,2),'') # set labels for x ticks
  1511. ticklabelsy <- ticklocsy # set labels for y ticks
  1513. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1514. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1515. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1516. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1517. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1518. title(xlab="Time (s)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1519. title(ylab="Diameter change (mm)", col.lab=rgb(0,0,0), line=2.8, cex.lab=1.5)
  1521. filt <- butter(3,0.083)
  1522. filtwaveP <- signal::filter(filt, meanwavesP)
  1523. filtU <- signal::filter(filt,meanwavesPCI[1,])
  1524. filtL <- signal::filter(filt,meanwavesPCI[2,])
  1525. polygon(times[c(1:length(times),length(times):1)],c(filtL,filtU[length(filtU):1]),col=addalpha(colpal[2],0.4),border=NA)
  1526. lines(times, filtwaveP, col=colpal[2], lwd=1.5)
  1527. lines(times, (sinewave/10)-0.9, col='black', lwd=2)
  1529. for (part in 1:24) {
  1530. load(file=paste0(d[part]))
  1531. filt <- butter(3,0.083)
  1532. filtwaveP <- signal::filter(filt, meanwavesP)
  1533. lines(times, filtwaveP, col=addalpha(colpal[2],0.7), lwd=0.2)
  1534. }
  1536. text(-1,0.4,'(p)',pos=4,cex=2)
  1538. load(file=paste0(datadir,'AveragedataLum.RData'))
  1540. plotlims <- c(0,2,0,0.15)
  1541. ticklocsx <- seq(0,2,0.5) # locations of tick marks on x axis
  1542. ticklocsy <- seq(0,0.15,0.05) # locations of tick marks on y axis
  1543. ticklabelsx <-ticklocsx # set labels for x ticks
  1544. ticklabelsy <- ticklocsy # set labels for y ticks
  1546. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1547. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1548. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1549. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1550. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1551. title(xlab="Frequency (Hz)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1552. title(ylab="Amplitude (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1554. polygon(frequencies[c(3:21,21:3)],c(meanspectraPCI[1,3:21],meanspectraPCI[2,21:3]),col=addalpha(colpal[2],0.3),border=NA)
  1556. lines(frequencies[3:21],abs(meanspectraP[3:21]), col=colpal[2], lwd=1.5, cex=0.5)
  1558. text(0,0.14,'(q)',pos=4,cex=2)
  1561. load(file=paste0(datadir,'AveragedataLM.RData'))
  1564. plotlims <- c(0,2,0,0.015)
  1565. ticklocsx <- seq(0,2,0.5) # locations of tick marks on x axis
  1566. ticklocsy <- seq(0,0.015,0.005) # locations of tick marks on y axis
  1567. ticklabelsx <-ticklocsx # set labels for x ticks
  1568. ticklabelsy <- ticklocsy # set labels for y ticks
  1570. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1571. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1572. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1573. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1574. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1575. title(xlab="Frequency (Hz)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1576. title(ylab="Amplitude (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1578. polygon(frequencies[c(3:21,21:3)],c(meanspectraPCI[1,3:21],meanspectraPCI[2,21:3]),col=addalpha(colpal[2],0.3),border=NA)
  1580. lines(frequencies[3:21],abs(meanspectraP[3:21]), col=colpal[2], lwd=1.5, cex=0.5)
  1582. text(0,0.014,'(r)',pos=4,cex=2)
  1585. load(file=paste0(datadir,'AveragedataS.RData'))
  1587. plotlims <- c(0,2,0,0.015)
  1588. ticklocsx <- seq(0,2,0.5) # locations of tick marks on x axis
  1589. ticklocsy <- seq(0,0.015,0.005) # locations of tick marks on y axis
  1590. ticklabelsx <-ticklocsx # set labels for x ticks
  1591. ticklabelsy <- ticklocsy # set labels for y ticks
  1593. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1594. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1595. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1596. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1597. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1598. title(xlab="Frequency (Hz)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1599. title(ylab="Amplitude (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1601. polygon(frequencies[c(3:21,21:3)],c(meanspectraPCI[1,3:21],meanspectraPCI[2,21:3]),col=addalpha(colpal[2],0.3),border=NA)
  1603. lines(frequencies[3:21],abs(meanspectraP[3:21]), col=colpal[2], lwd=1.5, cex=0.5)
  1605. text(0,0.014,'(s)',pos=4,cex=2)
  1608. load(file=paste0(datadir,'AveragedataMel.RData'))
  1610. plotlims <- c(0,2,0,0.015)
  1611. ticklocsx <- seq(0,2,0.5) # locations of tick marks on x axis
  1612. ticklocsy <- seq(0,0.015,0.005) # locations of tick marks on y axis
  1613. ticklabelsx <-ticklocsx # set labels for x ticks
  1614. ticklabelsy <- ticklocsy # set labels for y ticks
  1616. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1617. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1618. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1619. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1620. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1621. title(xlab="Frequency (Hz)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5)
  1622. title(ylab="Amplitude (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1624. polygon(frequencies[c(3:21,21:3)],c(meanspectraPCI[1,3:21],meanspectraPCI[2,21:3]),col=addalpha(colpal[2],0.3),border=NA)
  1626. lines(frequencies[3:21],abs(meanspectraP[3:21]), col=colpal[2], lwd=1.5, cex=0.5)
  1628. text(0,0.014,'(t)',pos=4,cex=2)
  1631. dev.off()
  1633. }
  1635. ```
  1637. ```{r plot1Fdata, include=FALSE, results='hide'}
  1639. if (processdata > 0){
  1641. contrastsfinedB <- 2:40
  1643. pdf(paste0(figdir,"FirstHarmonic.pdf"), bg="transparent", height = 24, width = 8)
  1645. layout(matrix(c(1,2,3,4,5,5,5,5,6,7,8,9),6,2,byrow=TRUE))
  1647. par(mar=c(3,6,2,1))
  1649. condlist <- c(1,2,6)
  1650. condlist2 <- c(1,2,4)
  1651. plotlims <- c(15,40,0,0.15)
  1652. plotlims2 <- c(15,40,0,0.015)
  1653. ticklocsx <- contrastsdB # locations of tick marks on x axis
  1654. ticklocsy <- seq(0,0.15,0.05)
  1655. ticklocsy2 <- seq(0,0.015,0.005) # locations of tick marks on y axis
  1656. ticklabelsx <-c(6,12,24,48,96) # set labels for x ticks
  1657. ticklabelsy <- ticklocsy
  1658. ticklabelsy2 <- ticklocsy2 # set labels for y ticks
  1660. load(file=paste0(datadir,'AveragedataLum.RData'))
  1662. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1663. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1664. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1665. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1666. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1667. title(xlab="Target contrast (%)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1668. title(ylab="Amplitude at 0.5Hz (mm)", col.lab=rgb(0,0,0), line=3, cex.lab=1.5)
  1669. title(main='Light flux',cex.main=2)
  1670. for (cond in 1:3){
  1671. lines(contrastsfinedB,modrespsLum[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1672. # lines(contrastsfinedB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1673. arrows(contrastsdB,abs(cleanmeansP[condlist[cond],]),contrastsdB,cleanmeansPCI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1674. arrows(contrastsdB,abs(cleanmeansP[condlist[cond],]),contrastsdB,cleanmeansPCI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1675. points(contrastsdB,abs(cleanmeansP[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1676. }
  1677. legend(15,0.125,c('Monocular','Binocular','Dichoptic cross'),pch=21:23,pt.bg=colpal[condlist],pt.lwd=3,pt.cex=1.6,box.lwd=2,cex=1.2)
  1678. text(17,0.14,'(a)',adj=0.5,cex=2.5)
  1681. load(file=paste0(datadir,'AveragedataLM.RData'))
  1683. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  1684. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1685. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  1686. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1687. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  1688. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1689. title(ylab="Amplitude at 0.5Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1690. title(main='L-M cone',cex.main=2)
  1691. for (cond in 1:3){
  1692. lines(contrastsfinedB,modrespsLM[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1693. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1694. arrows(contrastsdB,abs(cleanmeansP[condlist[cond],]),contrastsdB,cleanmeansPCI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1695. arrows(contrastsdB,abs(cleanmeansP[condlist[cond],]),contrastsdB,cleanmeansPCI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1696. points(contrastsdB,abs(cleanmeansP[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1697. }
  1698. text(17,0.014,'(b)',adj=0.5,cex=2.5)
  1701. load(file=paste0(datadir,'AveragedataS.RData'))
  1703. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  1704. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1705. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  1706. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1707. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  1708. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1709. title(ylab="Amplitude at 0.5Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1710. title(main='S cone',cex.main=2)
  1711. for (cond in 1:3){
  1712. lines(contrastsfinedB,modrespsS[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1713. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1714. arrows(contrastsdB,abs(cleanmeansP[condlist[cond],]),contrastsdB,cleanmeansPCI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1715. arrows(contrastsdB,abs(cleanmeansP[condlist[cond],]),contrastsdB,cleanmeansPCI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1716. points(contrastsdB,abs(cleanmeansP[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1717. }
  1718. text(17,0.014,'(c)',adj=0.5,cex=2.5)
  1721. load(file=paste0(datadir,'AveragedataMel.RData'))
  1723. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  1724. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1725. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  1726. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1727. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  1728. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1729. title(ylab="Amplitude at 0.5Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1730. title(main='Melanopsin',cex.main=2)
  1731. for (cond in 1:3){
  1732. lines(contrastsfinedB,modrespsMel[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1733. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1734. arrows(contrastsdB,abs(cleanmeansP[condlist[cond],]),contrastsdB,cleanmeansPCI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1735. arrows(contrastsdB,abs(cleanmeansP[condlist[cond],]),contrastsdB,cleanmeansPCI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1736. points(contrastsdB,abs(cleanmeansP[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1737. }
  1738. text(17,0.014,'(d)',adj=0.5,cex=2.5)
  1741. par(mar=c(8,12,8,4))
  1743. plotlims <- c(-24,12,0,0.6)
  1744. ticklocsx <- seq(-24,12,12)
  1745. ticklocsy <- seq(0,0.6,0.1)
  1746. ticklabelsx <- c('1/16','1/4','1','4')
  1747. ticklabelsy <- ticklocsy # set labels for y ticks
  1749. par(pty="s") # make axis square
  1750. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4]) # create an empty axis of the correct dimensions
  1751. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2)
  1752. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1753. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=1, cex=1.5)
  1754. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1, cex=1.5)
  1756. title(xlab=expression(omega), col.lab=rgb(0,0,0), line=4, cex.lab=6)
  1757. title(ylab="Probability", col.lab=rgb(0,0,0), line=3, cex.lab=3)
  1758. title(main='Suppressive weight 1F',cex.main=3)
  1760. a <- density(20*log10(pupilweightpopLum),from=-48,to=12)
  1761. polygon(c(-48,a$x,12), c(0,a$y,0), col=addalpha(stimcolours[1],0.5),border=NA)
  1763. a<- density(20*log10(pupilweightpopLM),from=-48,to=12)
  1764. polygon(c(-48,a$x,12), c(0,a$y,0), col=addalpha(stimcolours[4],0.5),border=NA)
  1766. a <- density(20*log10(pupilweightpopS),from=-48,to=12)
  1767. polygon(c(-48,a$x,12), c(0,a$y,0), col=addalpha(stimcolours[6],0.5),border=NA)
  1769. a <- density(20*log10(pupilweightpopMel),from=-48,to=12)
  1770. polygon(c(-48,a$x,12), c(0,a$y,0), col=addalpha(stimcolours[7],0.5),border=NA)
  1772. prior <- dnorm(10^((-48:12)/20),1,0.5)
  1773. prior <- 0.5*prior/max(prior)
  1774. lines(-48:12,prior,lwd=2)
  1776. lines(20*log10(c(mean(pupilWLum),mean(pupilWLum))),c(0,0.6),lwd=6,col=stimcolours[1])
  1777. lines(20*log10(c(mean(pupilWLM),mean(pupilWLM))),c(0,0.6),lwd=6,col=stimcolours[4])
  1778. lines(20*log10(c(mean(pupilWS),mean(pupilWS))),c(0,0.6),lwd=6,col=stimcolours[6])
  1779. lines(20*log10(c(mean(pupilWMel),mean(pupilWMel))),c(0,0.6),lwd=6,col=stimcolours[7])
  1781. legend(-24,0.6,c('Prior','LF 1F', 'L-M 1F', 'S 1F', 'Mel 1F'), col=c('black',stimcolours[c(1,4,6,7)]),lwd=c(2,6,6,6,6),lty=1, cex=1.8, box.lwd=2, bg='white')
  1783. text(10,0.58,'(e)',adj=0.5,cex=2.5)
  1786. par(mar=c(3,6,2,1))
  1788. condlist <- c(4,5,6)
  1789. condlist2 <- c(1,4,5)
  1790. plotlims <- c(15,40,0,0.15)
  1791. plotlims2 <- c(15,40,0,0.015)
  1792. ticklocsx <- contrastsdB # locations of tick marks on x axis
  1793. ticklocsy <- seq(0,0.15,0.05)
  1794. ticklocsy2 <- seq(0,0.015,0.005) # locations of tick marks on y axis
  1795. ticklabelsx <-c(6,12,24,48,96) # set labels for x ticks
  1796. ticklabelsy <- ticklocsy
  1797. ticklabelsy2 <- ticklocsy2 # set labels for y ticks
  1799. load(file=paste0(datadir,'AveragedataLum.RData'))
  1801. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4])
  1802. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1803. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  1804. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1805. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1)
  1806. title(xlab="Target contrast (%)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1807. title(ylab="Amplitude at 0.4Hz (mm)", col.lab=rgb(0,0,0), line=3, cex.lab=1.5)
  1808. title(main='Light flux',cex.main=2)
  1809. for (cond in 1:3){
  1810. lines(contrastsfinedB,modrespsLum[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1811. # lines(contrastsfinedB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1812. arrows(contrastsdB,abs(cleanmasksP[condlist[cond],]),contrastsdB,cleanmasksPCI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1813. arrows(contrastsdB,abs(cleanmasksP[condlist[cond],]),contrastsdB,cleanmasksPCI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1814. points(contrastsdB,abs(cleanmasksP[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1815. }
  1816. legend(15,0.125,c('Monocular','Binocular cross','Dichoptic cross'),pch=21:23,pt.bg=colpal[condlist],pt.lwd=3,pt.cex=1.6,box.lwd=2,cex=1.2)
  1817. text(17,0.14,'(f)',adj=0.5,cex=2.5)
  1820. load(file=paste0(datadir,'AveragedataLM.RData'))
  1822. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  1823. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1824. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  1825. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1826. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  1827. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1828. title(ylab="Amplitude at 0.4Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1829. title(main='L-M cone',cex.main=2)
  1830. for (cond in 1:3){
  1831. lines(contrastsfinedB,modrespsLM[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1832. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1833. arrows(contrastsdB,abs(cleanmasksP[condlist[cond],]),contrastsdB,cleanmasksPCI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1834. arrows(contrastsdB,abs(cleanmasksP[condlist[cond],]),contrastsdB,cleanmasksPCI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1835. points(contrastsdB,abs(cleanmasksP[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1836. }
  1837. text(17,0.014,'(g)',adj=0.5,cex=2.5)
  1840. load(file=paste0(datadir,'AveragedataS.RData'))
  1842. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  1843. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1844. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  1845. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1846. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  1847. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1848. title(ylab="Amplitude at 0.4Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1849. title(main='S cone',cex.main=2)
  1850. for (cond in 1:3){
  1851. lines(contrastsfinedB,modrespsS[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1852. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1853. arrows(contrastsdB,abs(cleanmasksP[condlist[cond],]),contrastsdB,cleanmasksPCI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1854. arrows(contrastsdB,abs(cleanmasksP[condlist[cond],]),contrastsdB,cleanmasksPCI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1855. points(contrastsdB,abs(cleanmasksP[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1856. }
  1857. text(17,0.014,'(h)',adj=0.5,cex=2.5)
  1860. load(file=paste0(datadir,'AveragedataMel.RData'))
  1862. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  1863. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1864. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  1865. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1866. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  1867. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1868. title(ylab="Amplitude at 0.4Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1869. title(main='Melanopsin',cex.main=2)
  1870. for (cond in 1:3){
  1871. lines(contrastsfinedB,modrespsMel[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1872. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1873. arrows(contrastsdB,abs(cleanmasksP[condlist[cond],]),contrastsdB,cleanmasksPCI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1874. arrows(contrastsdB,abs(cleanmasksP[condlist[cond],]),contrastsdB,cleanmasksPCI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1875. points(contrastsdB,abs(cleanmasksP[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1876. }
  1877. text(17,0.014,'(i)',adj=0.5,cex=2.5)
  1881. dev.off()
  1883. }
  1885. ```
  1887. ```{r plot2Fdata, include=FALSE, results='hide'}
  1889. if (processdata > 0){
  1891. contrastsfinedB <- 2:40
  1893. pdf(paste0(figdir,"SecondHarmonic.pdf"), bg="transparent", height = 24, width = 8)
  1895. layout(matrix(c(1,2,3,4,5,5,5,5,6,7,8,9),6,2,byrow=TRUE))
  1896. par(mar=c(3,6,2,1))
  1898. condlist <- c(1,2,6)
  1899. condlist2 <- c(1,2,4)
  1900. plotlims <- c(15,40,0,0.15)
  1901. plotlims2 <- c(15,40,0,0.015)
  1902. ticklocsx <- contrastsdB # locations of tick marks on x axis
  1903. ticklocsy <- seq(0,0.15,0.05)
  1904. ticklocsy2 <- seq(0,0.015,0.005) # locations of tick marks on y axis
  1905. ticklabelsx <-c(6,12,24,48,96) # set labels for x ticks
  1906. ticklabelsy <- ticklocsy
  1907. ticklabelsy2 <- ticklocsy2 # set labels for y ticks
  1909. load(file=paste0(datadir,'AveragedataLum.RData'))
  1911. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  1912. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1913. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  1914. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1915. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  1916. title(xlab="Target contrast (%)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1917. title(ylab="Amplitude at 1Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1918. title(main='Light flux',cex.main=2)
  1919. for (cond in 1:3){
  1920. lines(contrastsfinedB,modrespsLum2[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1921. # lines(contrastsfinedB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1922. arrows(contrastsdB,abs(cleanmeansP2[condlist[cond],]),contrastsdB,cleanmeansP2CI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1923. arrows(contrastsdB,abs(cleanmeansP2[condlist[cond],]),contrastsdB,cleanmeansP2CI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1924. points(contrastsdB,abs(cleanmeansP2[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1925. }
  1926. legend(15,0.0125,c('Monocular','Binocular','Dichoptic cross'),pch=21:23,pt.bg=colpal[condlist],pt.lwd=3,pt.cex=1.6,box.lwd=2,cex=1.2)
  1927. text(17,0.014,'(a)',adj=0.5,cex=2.5)
  1930. load(file=paste0(datadir,'AveragedataLM.RData'))
  1932. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  1933. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1934. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  1935. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1936. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  1937. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1938. title(ylab="Amplitude at 1Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1939. title(main='L-M cone',cex.main=2)
  1940. for (cond in 1:3){
  1941. lines(contrastsfinedB,modrespsLM2[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1942. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1943. arrows(contrastsdB,abs(cleanmeansP2[condlist[cond],]),contrastsdB,cleanmeansP2CI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1944. arrows(contrastsdB,abs(cleanmeansP2[condlist[cond],]),contrastsdB,cleanmeansP2CI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1945. points(contrastsdB,abs(cleanmeansP2[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1946. }
  1947. text(17,0.014,'(b)',adj=0.5,cex=2.5)
  1950. load(file=paste0(datadir,'AveragedataS.RData'))
  1952. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  1953. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1954. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  1955. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1956. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  1957. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1958. title(ylab="Amplitude at 1Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1959. title(main='S cone',cex.main=2)
  1960. for (cond in 1:3){
  1961. lines(contrastsfinedB,modrespsS2[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1962. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1963. arrows(contrastsdB,abs(cleanmeansP2[condlist[cond],]),contrastsdB,cleanmeansP2CI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1964. arrows(contrastsdB,abs(cleanmeansP2[condlist[cond],]),contrastsdB,cleanmeansP2CI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1965. points(contrastsdB,abs(cleanmeansP2[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1966. }
  1967. text(17,0.014,'(c)',adj=0.5,cex=2.5)
  1970. load(file=paste0(datadir,'AveragedataMel.RData'))
  1972. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  1973. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  1974. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  1975. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  1976. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  1977. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  1978. title(ylab="Amplitude at 1Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  1979. title(main='Melanopsin',cex.main=2)
  1980. for (cond in 1:3){
  1981. lines(contrastsfinedB,modrespsMel2[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1982. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  1983. arrows(contrastsdB,abs(cleanmeansP2[condlist[cond],]),contrastsdB,cleanmeansP2CI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1984. arrows(contrastsdB,abs(cleanmeansP2[condlist[cond],]),contrastsdB,cleanmeansP2CI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  1985. points(contrastsdB,abs(cleanmeansP2[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  1986. }
  1987. text(17,0.014,'(d)',adj=0.5,cex=2.5)
  1990. par(mar=c(8,12,8,4))
  1992. plotlims <- c(-24,12,0,0.6)
  1993. ticklocsx <- seq(-24,12,12)
  1994. ticklocsy <- seq(0,0.6,0.1)
  1995. ticklabelsx <- c('1/16','1/4','1','4')
  1996. ticklabelsy <- ticklocsy # set labels for y ticks
  1998. par(pty="s") # make axis square
  1999. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims[1:2], ylim=plotlims[3:4]) # create an empty axis of the correct dimensions
  2000. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2)
  2001. axis(2, at=ticklocsy, tck=0.01, lab=F, lwd=2)
  2002. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=1, cex=1.5)
  2003. mtext(text = ticklabelsy, side = 2, at=ticklocsy, line=0.2, las=1, cex=1.5)
  2005. title(xlab=expression(omega), col.lab=rgb(0,0,0), line=4, cex.lab=6)
  2006. title(ylab="Probability", col.lab=rgb(0,0,0), line=3, cex.lab=3)
  2007. title(main='Suppressive weight 2F',cex.main=3)
  2009. a <- density(20*log10(pupilweightpopLum2),from=-48,to=12)
  2010. polygon(c(-48,a$x,12), c(0,a$y,0), col=addalpha(stimcolours[1],0.5),border=NA)
  2012. a<- density(20*log10(pupilweightpopLM2),from=-48,to=12)
  2013. polygon(c(-48,a$x,12), c(0,a$y,0), col=addalpha(stimcolours[4],0.5),border=NA)
  2015. a <- density(20*log10(pupilweightpopS2),from=-48,to=12)
  2016. polygon(c(-48,a$x,12), c(0,a$y,0), col=addalpha(stimcolours[6],0.5),border=NA)
  2018. a <- density(20*log10(pupilweightpopMel2),from=-48,to=12)
  2019. polygon(c(-48,a$x,12), c(0,a$y,0), col=addalpha(stimcolours[7],0.5),border=NA)
  2021. prior <- dnorm(10^((-48:12)/20),1,0.5)
  2022. prior <- 0.5*prior/max(prior)
  2023. lines(-48:12,prior,lwd=2)
  2025. lines(20*log10(c(mean(pupilWLum2),mean(pupilWLum2))),c(0,0.6),lwd=6,col=stimcolours[1])
  2026. lines(20*log10(c(mean(pupilWLM2),mean(pupilWLM2))),c(0,0.6),lwd=6,col=stimcolours[4])
  2027. lines(20*log10(c(mean(pupilWS2),mean(pupilWS2))),c(0,0.6),lwd=6,col=stimcolours[6])
  2028. lines(20*log10(c(mean(pupilWMel2),mean(pupilWMel2))),c(0,0.6),lwd=6,col=stimcolours[7])
  2030. legend(-24,0.6,c('Prior','LF 2F', 'L-M 2F', 'S 2F', 'Mel 2F'), col=c('black',stimcolours[c(1,4,6,7)]),lwd=c(2,6,6,6,6),lty=1, cex=1.8, box.lwd=2, bg='white')
  2032. text(10,0.58,'(e)',adj=0.5,cex=2.5)
  2035. par(mar=c(3,6,2,1))
  2037. condlist <- c(4,5,6)
  2038. condlist2 <- c(1,4,5)
  2039. plotlims <- c(15,40,0,0.15)
  2040. plotlims2 <- c(15,40,0,0.015)
  2041. ticklocsx <- contrastsdB # locations of tick marks on x axis
  2042. ticklocsy <- seq(0,0.15,0.05)
  2043. ticklocsy2 <- seq(0,0.015,0.005) # locations of tick marks on y axis
  2044. ticklabelsx <-c(6,12,24,48,96) # set labels for x ticks
  2045. ticklabelsy <- ticklocsy
  2046. ticklabelsy2 <- ticklocsy2 # set labels for y ticks
  2048. load(file=paste0(datadir,'AveragedataLum.RData'))
  2050. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  2051. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  2052. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  2053. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  2054. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  2055. title(xlab="Target contrast (%)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  2056. title(ylab="Amplitude at 0.8Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  2057. title(main='Light flux',cex.main=2)
  2058. for (cond in 1:3){
  2059. lines(contrastsfinedB,modrespsLum2[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  2060. # lines(contrastsfinedB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  2061. arrows(contrastsdB,abs(cleanmasksP2[condlist[cond],]),contrastsdB,cleanmasksP2CI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  2062. arrows(contrastsdB,abs(cleanmasksP2[condlist[cond],]),contrastsdB,cleanmasksP2CI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  2063. points(contrastsdB,abs(cleanmasksP2[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  2064. }
  2065. legend(15,0.0125,c('Monocular','Binocular cross','Dichoptic cross'),pch=21:23,pt.bg=colpal[condlist],pt.lwd=3,pt.cex=1.6,box.lwd=2,cex=1.2)
  2066. text(17,0.014,'(f)',adj=0.5,cex=2.5)
  2069. load(file=paste0(datadir,'AveragedataLM.RData'))
  2071. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  2072. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  2073. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  2074. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  2075. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  2076. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  2077. title(ylab="Amplitude at 0.8Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  2078. title(main='L-M cone',cex.main=2)
  2079. for (cond in 1:3){
  2080. lines(contrastsfinedB,modrespsLM2[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  2081. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  2082. arrows(contrastsdB,abs(cleanmasksP2[condlist[cond],]),contrastsdB,cleanmasksP2CI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  2083. arrows(contrastsdB,abs(cleanmasksP2[condlist[cond],]),contrastsdB,cleanmasksP2CI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  2084. points(contrastsdB,abs(cleanmasksP2[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  2085. }
  2086. text(17,0.014,'(g)',adj=0.5,cex=2.5)
  2089. load(file=paste0(datadir,'AveragedataS.RData'))
  2091. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  2092. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  2093. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  2094. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  2095. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  2096. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  2097. title(ylab="Amplitude at 0.8Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  2098. title(main='S cone',cex.main=2)
  2099. for (cond in 1:3){
  2100. lines(contrastsfinedB,modrespsS2[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  2101. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  2102. arrows(contrastsdB,abs(cleanmasksP2[condlist[cond],]),contrastsdB,cleanmasksP2CI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  2103. arrows(contrastsdB,abs(cleanmasksP2[condlist[cond],]),contrastsdB,cleanmasksP2CI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  2104. points(contrastsdB,abs(cleanmasksP2[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  2105. }
  2106. text(17,0.014,'(h)',adj=0.5,cex=2.5)
  2109. load(file=paste0(datadir,'AveragedataMel.RData'))
  2111. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=plotlims2[1:2], ylim=plotlims2[3:4])
  2112. axis(1, at=ticklocsx, tck=0.01, lab=F, lwd=2) # plot tick marks (no labels)
  2113. axis(2, at=ticklocsy2, tck=0.01, lab=F, lwd=2)
  2114. mtext(text = ticklabelsx, side = 1, at=ticklocsx, line=0.3) # add the tick labels
  2115. mtext(text = ticklabelsy2, side = 2, at=ticklocsy2, line=0.2, las=1)
  2116. title(xlab="Target contrast (% of max)", col.lab=rgb(0,0,0), line=1.5, cex.lab=1.5) # titles for axes
  2117. title(ylab="Amplitude at 0.8Hz (mm)", col.lab=rgb(0,0,0), line=3.5, cex.lab=1.5)
  2118. title(main='Melanopsin',cex.main=2)
  2119. for (cond in 1:3){
  2120. lines(contrastsfinedB,modrespsMel2[condlist2[cond],], col=colpal[condlist[cond]], lwd=3, cex=0.5)
  2121. # lines(contrastsdB,abs(cleanmeansP[condlist[cond],]), col=colpal[condlist[cond]], lwd=3, cex=0.5)
  2122. arrows(contrastsdB,abs(cleanmasksP2[condlist[cond],]),contrastsdB,cleanmasksP2CI[1,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  2123. arrows(contrastsdB,abs(cleanmasksP2[condlist[cond],]),contrastsdB,cleanmasksP2CI[2,condlist[cond],],lwd=3,col=colpal[condlist[cond]],length=0.05,angle=90)
  2124. points(contrastsdB,abs(cleanmasksP2[condlist[cond],]), pch = 20+cond, col='black', bg=colpal[condlist[cond]], cex=1.6, lwd=3)
  2125. }
  2126. text(17,0.014,'(i)',adj=0.5,cex=2.5)
  2130. dev.off()
  2132. }
  2134. ```
  2136. ```{r makephaseplots, include=FALSE, results='hide'}
  2138. # generate plots showing the phase angles
  2139. if (processdata > 0){
  2141. load(file=paste0(datadir,'AveragedataLum.RData'))
  2142. xlum <- Re(cleanmeansP)
  2143. ylum <- Im(cleanmeansP)
  2144. xlum2 <- Re(cleanmeansP2)
  2145. ylum2 <- Im(cleanmeansP2)
  2146. load(file=paste0(datadir,'AveragedataMel.RData'))
  2147. xmel <- Re(cleanmeansP)
  2148. ymel <- Im(cleanmeansP)
  2149. xmel2 <- Re(cleanmeansP2)
  2150. ymel2 <- Im(cleanmeansP2)
  2151. load(file=paste0(datadir,'AveragedataLM.RData'))
  2152. xlm <- Re(cleanmeansP)
  2153. ylm <- Im(cleanmeansP)
  2154. xlm2 <- Re(cleanmeansP2)
  2155. ylm2 <- Im(cleanmeansP2)
  2156. load(file=paste0(datadir,'AveragedataS.RData'))
  2157. xs <- Re(cleanmeansP)
  2158. ys <- Im(cleanmeansP)
  2159. xs2 <- Re(cleanmeansP2)
  2160. ys2 <- Im(cleanmeansP2)
  2162. pdf(paste(figdir,"phaseplots.pdf",sep=''), bg="transparent", height = 5, width = 10)
  2164. par(mfrow=c(1,2),mar=c(1,1,1,1))
  2167. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=c(-0.022,0.022), ylim=c(-0.022,0.022))
  2169. rangles <- seq(0,360)*(pi/180)
  2170. r <- seq(0.005,0.02,0.005)
  2171. for (n in 1:length(r)){
  2172. lines(r[n]*cos(rangles),r[n]*sin(rangles),col='gray')
  2173. }
  2174. lines(c(0,0),c(-0.02,0.02),lwd=2)
  2175. lines(c(-0.02,0.02),c(0,0),lwd=2)
  2177. text(0.015,-0.001,'Amplitude',cex=1.25,adj=0.5)
  2178. text(0.021,0,'Phase',cex=1.25,adj=0.5,srt=-90)
  2180. lines(xlum[2,]/10,ylum[2,]/10,col=stimcolours[1],lwd=3)
  2181. points(xlum[2,]/10,ylum[2,]/10,pch=21,col=stimcolours[1],bg=stimcolours[2],cex=2,lwd=3)
  2183. lines(xmel[2,],ymel[2,],col=stimcolours[7],lwd=3)
  2184. points(xmel[2,],ymel[2,],pch=21,col=stimcolours[7],bg=stimcolours[8],cex=2,lwd=3)
  2186. lines(xlm[2,],ylm[2,],col=stimcolours[3],lwd=3)
  2187. points(xlm[2,],ylm[2,],pch=21,col=stimcolours[3],bg=stimcolours[4],cex=2,lwd=3)
  2189. lines(xs[2,],ys[2,],col=stimcolours[6],lwd=3)
  2190. points(xs[2,],ys[2,],pch=21,col=stimcolours[6],bg=stimcolours[5],cex=2,lwd=3)
  2192. text(-0.019,0.019,'(a)',cex=3)
  2193. text(0.016,0.019,'1F',cex=2.5)
  2195. legend(-0.019,0.0155,c('Light flux','L-M','S','Melanopsin'),lty=1,col=stimcolours[c(1,3,6,7)],lwd=3,pch=21,pt.cex=2,pt.bg=stimcolours[c(2,4,5,8)],cex=1.1,box.lwd=3,bg='white')
  2198. plot(x=NULL,y=NULL,axes=FALSE, ann=FALSE, xlim=c(-0.022,0.022), ylim=c(-0.022,0.022))
  2200. rangles <- seq(0,360)*(pi/180)
  2201. r <- seq(0.005,0.02,0.005)
  2202. for (n in 1:length(r)){
  2203. lines(r[n]*cos(rangles),r[n]*sin(rangles),col='gray')
  2204. }
  2205. lines(c(0,0),c(-0.02,0.02),lwd=2)
  2206. lines(c(-0.02,0.02),c(0,0),lwd=2)
  2208. text(0.015,-0.001,'Amplitude',cex=1.25,adj=0.5)
  2209. text(0.021,0,'Phase',cex=1.25,adj=0.5,srt=-90)
  2211. lines(xlum2[2,],ylum2[2,],col=stimcolours[1],lwd=3)
  2212. points(xlum2[2,],ylum2[2,],pch=21,col=stimcolours[1],bg=stimcolours[2],cex=2,lwd=3)
  2214. lines(xmel2[2,],ymel2[2,],col=stimcolours[7],lwd=3)
  2215. points(xmel2[2,],ymel2[2,],pch=21,col=stimcolours[7],bg=stimcolours[8],cex=2,lwd=3)
  2217. lines(xlm2[2,],ylm2[2,],col=stimcolours[3],lwd=3)
  2218. points(xlm2[2,],ylm2[2,],pch=21,col=stimcolours[3],bg=stimcolours[4],cex=2,lwd=3)
  2220. lines(xs2[2,],ys2[2,],col=stimcolours[6],lwd=3)
  2221. points(xs2[2,],ys2[2,],pch=21,col=stimcolours[6],bg=stimcolours[5],cex=2,lwd=3)
  2224. text(-0.019,0.019,'(b)',cex=3)
  2225. text(0.016,0.019,'2F',cex=2.5)
  2229. dev.off()
  2231. }
  2233. ```
  2235. # Abstract
  2237. Pupil diameters are regulated by the autonomic nervous system, which combines light signals across the eyes independently of the visual cortex. Distinct classes of retinal photoreceptor are involved in this process, with cones and rods driving the initial constriction and intrinsically photosensitive retinal ganglion cells maintaining diameter over prolonged time periods. We investigated binocular combination by targeting different photoreceptor pathways using a novel binocular multiprimary system to modulate the input spectra via silent substitution. At the first harmonic of the modulation frequency, luminance and S-cone responses showed strong binocular facilitation, and weak interocular suppression. Melanopsin responses were invariant to the number of eyes stimulated. The L-M pathway involved binocular inhibition, whereby responses to binocular stimulation were weaker than for monocular stimulation. The second harmonic involved strong interocular suppression in all pathways, but with some evidence of binocular facilitation. Our results are consistent with a computational model of binocular signal combination (implemented in a Bayesian hierarchical framework), in which the weight of interocular suppression differs across pathways. We also find pathway differences in response phase, consistent with different lag times for phototransduction. These results provide evidence that binocular interactions in the pupillary pathway differ across photoreceptor-directed modulations and across harmonic components of the response.
  2239. # Introduction
  2241. The autonomic nervous system regulates many involuntary bodily processes, including the constriction and dilation of the pupils in response to light [@McDougal2015]. The anatomical pathway from the retina to the subcortical nuclei controlling the pupillary light response (PLR) is well established: it includes the olivary pretectal nucleus (OPN), the Superior Cervical ganglion and the Edinger-Westphal nucleus, which project to the iris sphincter muscles that directly control the pupil size [@McDougal2015]. Evidence of a binocular component to the PLR is shown by the consensual response of the pupil (stimulation of one eye will cause constriction of the other eye) [@Wyatt1981]. The anatomical segregation of the subcortical pathway from the rest of the brain means that this binocular combination of signals must occur independently of the cortical processes of binocular integration required for visual perception. Our recent work [@Segala2023] has shown that the algorithm underlying binocular combination of light in the pupil pathway differs from that in the cortex. Here we extend this paradigm to compare binocular combination across signals from different photoreceptor classes that feed into the autonomic nervous system.
  2243. Different classes of retinal photoreceptors, including cones, rods and melanopsin-containing intrinsically photosensitive retinal ganglion cells (ipRGCs), are directly involved in controlling and maintaining the size of the pupils [e.g. @McDougal2010; @Barrionuevo2016; @Murray2018; @Spitschan2014; @Woelders2018]. Cones drive the initial rapid constriction of the pupils [@Mathot2018], while the slower and longer activation of the ipRGCs maintains constriction over a prolonged period of time and regulates the post-illumination pupillary response [@Markwell2010; @McDougal2010]. The ipRGCs are a recently discovered photoreceptor class [@Provencio2000] that express the photopigment melanopsin, and are involved in the regulation of the circadian rhythm [@Panda2002; @Ruby2002], forming a major input to the OPN [@Dacey2003]. The first direct evidence of the involvement of the ipRGCs in the PLR was shown in melanopsin knockout mice [@Lucas2003], resulting in the loss of the intrinsic photosensitivity of the cells and a reduced pupil constriction. Similar behaviour was later observed in primates and humans [@Gamlin2007] using silent substitution (see Methods), where it was shown that the PLR continues during light presentation even when cone and rod signalling is blocked, indicating the primary role of the ipRGCs in maintaining pupil constriction over a prolonged time.
  2245. Binocular combination has been extensively studied in visual perception, where it is mediated by neurons in primary visual cortex [@Hubel1962]. For pattern vision, binocular summation occurs at threshold, such that a lower contrast is required to detect a target shown to both eyes than a target shown to one eye [@Baker2018; @Campbell1965]. At higher contrasts, the phenomenon of 'ocularity invariance' is observed, in which the response to monocularly- and binocularly-presented patterns is equal [@Baker2007; @Moradi2009]. This is explained by a process of interocular suppression that cancels out the additional excitatory drive caused by stimulating two eyes. Our recent work [@Segala2023] showed that cortical signal combination for luminance flicker is substantially more linear than for spatial patterns, whereas there is evidence of interocular suppression in the pupil pathway. Here we measure the amplitude of pupil modulations in response to flickering stimuli presented as light flux, or directed towards the L-M cones, S-cones and melanopsin (see Figure \@ref(fig:spectraplots)e-l). Our design uses three comparisons, each estimating a distinct quantity. To obtain an objective measure of pupil responses, we presented the two eyes with slightly different temporal frequencies (0.5 and 0.4 Hz). Responses at these two frequencies fall in separate bins of the Fourier spectrum and so this 'frequency tagging' lets us recover the response driven by one eye while the other eye supplies a concurrent, separable signal. Comparing monocular and binocular stimulation measures the net effect of adding a second eye, but cannot by itself separate additional excitation from interocular suppression. Therefore we used a dichoptic mask, in which the target eye is stimulated while the other eye views a fixed-contrast modulation. This allows us to isolate the suppression that one eye exerts on the response to the other. Together these comparisons allow us to test whether binocular facilitation and interocular suppression differ across photoreceptor-directed modulations. We interpret the results using a contemporary model of binocular vision [@Meese2006], fitted within a hierarchical Bayesian framework.
  2247. # Materials and methods
  2249. ## Participants
  2251. Twenty-four participants were recruited for each of the four experiments for a total of ninety-six adult participants (28 male, 68 female), whose ages ranged from 18 to 41. All participants had normal or corrected to normal vision, no known abnormalities of binocular or colour vision, and gave written informed consent. Our procedures were approved by the Ethics Committee of the Department of Psychology at the University of York (identification number 184). Participant recruitment began on 3 March 2023 and concluded on 8 Decemeber 2023.
  2253. ## Apparatus & stimuli
  2255. To present synchronised stimulus modulations independently to each eye, two light engines (SpectraTuneLAB, which are thermally stable across time with active cooling according to the manufacturers: LEDMOTIVE Technologies, LLC, Barcelona, Spain), each with 10 independently addressable LED colour channels, were integrated into a customised binocular viewing system. The light engines were operated via a Python interface to their REST API [@Martin2022], which supports synchronous launch and playback of spectral sequences prepared in advance and stored in JSON format. A special command was commissioned from LEDMOTIVE to allow two different stimulus files to be launched simultaneously from the two devices. The synchronisation of the spectra from the two devices was tested and showed that they were synchronised to within about 3ms.
  2257. When preparing the spectral sequences, the age of participants was used to account for the yellowing of the lenses. We used the silent substitution technique [@Estevez1982] to selectively stimulate specific photoreceptor classes. Silent substitution exploits the fact that each photoreceptor class has a distinct spectral tuning that overlaps with the others. Using a multiprimary system, in which the primaries (i.e. LEDs) have different spectra, it is possible to target one class of photoreceptors while maintaining the others at a constant activity level, effectively silencing them [@Shapiro1996; @Spitschan2018, this paper also offers a clear explanation of how to implement silent substitution]. We calculated silent substitution solutions using the _PySilSub_ toolbox [@Martin2023], using linear algebra. The outputs of the two light engines (see Figure \@ref(fig:spectraplots)b,c) were calibrated using an Ocean Optics Jaz spectroradiometer, which was wavelength-calibrated to an Argon lamp and intensity calibrated using a NIST-traceable light source. Each primary was also linearised using a polynomial fit (see Figure \@ref(fig:primarylinearity) for details). We used the 10-degree cone fundamentals [@Stockman2000], and estimates of melanopsin absorbance spectra from CIE S 026 (discussed in a previous paper [@Martin2023]) to calculate $\alpha$-opic irradiance.
  2259. ```{r spectraplots, fig.cap="Summary of the spectral power distributions and alpha-opic irradiances for the background and each condition, as well as averaged pupil diameters and Fourier spectra. Panel (a) shows a schematic of the binocular stimulation system for presenting spectrally tuned modulations independently to each eye. The VR headset was attached to a clamp stand that the experimenters could use to adjust the height and align the headset with the eyes of the participant. The participant's head was supported by a chin rest to keep it in position throughout the experiment. Panel (b) shows the outputs of each LED primary at maximum intensity, and panels (c) and (d) show the overall spectral power distributions and the alpha-opic irradiances of the background spectra used for both eyes. The subsequent rows show the power differences (e-h), and photoreceptor contrasts (i-l) relative to the background, averaged pupil diameter waveforms (m-p) and Fourier spectra (q-t) for binocular stimulation. Column headings indicate the pathway stimulated, and shaded regions in panels m-t indicate bootstrapped 95\\% confidence intervals.", fig.align="center", out.width = '80%', echo=FALSE}
  2261. knitr::include_graphics(paste0(figdir,'methodsfigrods.pdf'))
  2263. ```
  2265. The output from the light engines was directed through liquid light guides (LLG3-8H: Thorlabs Ltd, Cambridgeshire, UK) and diffused onto semi-opaque and highly diffusive white glass discs with a diameter of 50 mm for even illumination (34-473: Edmund Optics, York, UK). The light guide gaskets were butt-coupled to the light engine diffusers with threaded adapters (SM1A9, AD3LLG: Thorlabs Ltd, Cambridgeshire, UK) and the exiting ends of the light guides were mated with 51 mm depth optical cylinders (SM2L20: Thorlabs Ltd, Cambridgeshire, UK) via appropriately threaded adapters (AD3LLG, SM2A6: Thorlabs Ltd, Cambridgeshire, UK). The stimulus diffuser discs were retained at the front end of the optical cylinders about 51 mm from the light source, at which distance the output beam was sufficiently dispersed to afford even illumination of the diffuser when viewed from the front. To guarantee safe illumination levels, a circular neutral-density filter with the same diameter as the white glass discs (50 mm) and an optical density of 0.6 log units was placed in the optical path between the light source and the diffusers. A small circular piece of blackout material with a diameter of about 8 degrees (10 mm) was positioned centrally on the front of each diffuser disc to aid as a fusion lock, as a fixation point, and to occlude the fovea.
  2267. The diffuser discs were positioned in the objective planes of the lenses of a modified VR headset (SHINECON SC-G01, Dongguan Shinecon Industrial Co. Ltd., Guangdong, China), which was used by the participants to view the stimuli. The stimuli were two discs of flickering light with a diameter of about 30 degrees, which were fused together into a cyclopean percept resembling a donut-shaped ring of light, similar to that used in other studies [e.g., @Barrionuevo2016; @Murray2018; @Spitschan2014; @Spitschan2015; @Zele2018]. The VR headset modifications allowed for small adjustments to account for individual differences in interpupillary distance and focal length. The use of this set up allowed us to modulate the stimuli in three different ocular configurations, similar to the ones we used in our previous study [@Segala2023]: monocular, binocular and dichoptic. In the monocular configuration, the unstimulated eye still saw a non-flickering disc of mean light flux. A schematic of the stimulation system is shown in Figure \@ref(fig:spectraplots)a. Pupillometry data were collected using a binocular Pupil Core eye-tracker headset (Pupil Labs GmbH, Berlin, Germany [@Kassner2014]) running at 120 Hz, and the signals were recorded with the Pupil Capture software.
  2269. Our previous study [@Segala2023] used a temporal frequency of 2Hz for foveal luminance flicker, and recorded EEG data simultaneously with pupillometry. Initial pilot experiments indicated that this frequency was too high to elicit measurable responses when stimulating individual photoreceptor pathways. For all experiments, we therefore used a primary flicker frequency of 0.5Hz, as previous literature showed that this was slow enough to elicit a pupil response from all photoreceptor classes [@Spitschan2014]. We also focussed on only recording pupillometry data as this frequency would be too slow to elicit steady-state EEG responses [@Norcia2015].
  2271. For all experiments, sinusoidal temporal modulations were presented against the same background spectrum (matched between the eyes), which was used to achieve silent substitution in the three photoreceptor modulation experiments. The background spectra were defined by setting all channels to half maximum output for the brighter of the two devices (STLab 1, left eye) and then using the STLab 1/STLab 2 calibration ratio to find the equivalent settings for the companion device (STLab 2, right eye). The background spectrum illuminance was about 74 lux, or 68.5 cd/m$^2$. The spectral power distributions and $\alpha$-opic irradiances of the background spectra for both eyes are shown in Figure \@ref(fig:spectraplots)c-d.
  2273. Silent substitution stimuli were prepared and calibrated for each participant with custom Python software [@Martin2023] and Python scripts. Estimates of photoreceptor spectral sensitivities for each participant were constructed from the known photopigment absorbance spectra [@Stockman2000], taking account of the peak axial density of the respective photopigments, as well as lens [@Pokorny1987; @Stockman1999; @Stockman2000] and macular pigment density [@Bone1988; @Stockman1999], in accordance with the field size and age-dependent CIEPO06 observer model [@CIE2006]. The melanopic and rhodopic action spectra of the 32-year-old standard observer were taken from CIE S 026 [@CIE2018] and then adjusted for age-related lens transmittance with a spectral correction function, in line with the standard. Macular pigment correction was not applied to the rhodopic and melanopic action spectra because rods are not present at the fovea and ipRGCs sit above the retinal pigment layer [@Trieschmann2007].
  2275. In the light flux experiment, the stimulus intensity was increased and decreased relative to the background, which we expected to modulate all photoreceptor classes (see Figure \@ref(fig:spectraplots)e,i). In the L-M cone modulation experiment, we used silent substitution to increase the L-cone activity, and simultaneously decrease the M-cone activity, during the first half-cycle of the sine wave. In the second half-cycle the polarity of the modulation reversed (see Figure \@ref(fig:spectraplots)f,j). The maximum available L-M contrast was about 10\%. In the S-cone modulation experiment, we increased and decreased S-cone-directed signals, whilst keeping activity in the other photoreceptors constant (see Figure \@ref(fig:spectraplots)g,k). Our system allowed a maximum contrast of 45\%. Finally, in the melanopsin experiment, we modulated the activity of the melanopsin-containing intrinsically photoreceptive retinal ganglion cells, whilst keeping cone activity constant (see Figure \@ref(fig:spectraplots)h,l). The maximum available melanopsin contrast was 22\%. We assume that the activity of rods was constant at the high background luminance intensity used here, and so did not attempt to silence rod activity in any condition, as this would have greatly reduced the available dynamic range. Splatter on nominally silenced photoreceptors was very small (see Figure \@ref(fig:spectraplots)j-l), well below the levels that would be expected to generate measurable pupil modulations, although we can observe that the rods may not be saturated in the L-M and melanopsin-directed stimuli (Figure \@ref(fig:spectraplots)j and \@ref(fig:spectraplots)l) and could intrude in these two experiments. We also estimated activation of penumbral L and M cones [@Spitschan2015; @Barrionuevo2016; @Zele2019], which was minimal ($\le 1.5%$\% contrast; less than splatter on the open-field cones) for melanopsin-directed stimuli. We note that the temporal frequency of our modulation (0.5Hz) is well below the range where penumbral cone activation can elicit visible percepts, and that such percepts fade after around 1 second [@Spitschan2015], and do not typically affect pupil responses [@Spitschan2014].
  2277. ## Procedure
  2279. Before the start of each experiment, participants adjusted the objective planes of the lenses with the help of the experimenter until the stimulus was in focus and they perceived the two pieces of blackout material as one fused disc. Pupil responses to binocular temporal contrast modulations were examined in a factorial design that combined six ocular conditions and five temporal contrast levels: 6, 12, 24, 48 and 96% of the available dynamic range. This design, similar to that used in our previous studies [@Segala2023; @Baker2017; @Baker2020], was applied in four separate experiments, each with a different mode of photoreceptor stimulation. In the first three conditions, the discs flickered at 0.5 Hz, in either a monocular, binocular or dichoptic arrangement. In the dichoptic condition the non-target eye saw a flickering fixed contrast of 48% of the available dynamic range. In the remaining three conditions (the cross-frequency conditions) one eye’s disc flickered at 0.4 Hz, and the other eye’s disc flickered at 0.5 Hz. This included monocular responses at 0.4 Hz, as well as binocular (one eye sees each frequency at the target contrast) and dichoptic (target stimulus flickering at 0.5 Hz, mask contrast of 48% at 0.4 Hz in the other eye) arrangements. We counterbalanced presentation of the target stimulus across the left and right eyes. The full set of conditions is summarised in Table \@ref(tab:diagnosticstable).
  2281. ```{r diagnosticstable, echo=FALSE}
  2283. diagnostics <- data.frame(
  2284. Fit = c("Light flux 1F", "L-M cone 1F", "S cone 1F", "Melanopsin 1F",
  2285. "Light flux 2F", "L-M cone 2F", "S cone 2F", "Melanopsin 2F"),
  2286. Rhat = c(1.00, 1.01, 1.01, 1.06, 1.01, 1.02, 1.02, 1.00),
  2287. ESS = c(17700, 9800, 6200, 10300, 5700, 3100, 3000, 11300),
  2288. RMSE = c(0.0043, 0.0016, 0.0028, 0.0030, 0.0014, 0.0021, 0.0020, 0.0035)
  2289. )
  2291. kable(diagnostics,
  2292. caption = "Convergence diagnostics and fit statistics for the eight model fits. Maximum split-$\\hat{R}$ is taken across all parameters; minimum effective sample size (ESS) is taken across the group-level parameters; RMSE is the root-mean-square error between the group-mean data and the model prediction.",
  2293. col.names = c("Fit", "max $\\hat{R}$", "min ESS", "RMSE (mm)"),
  2294. align = "lrrr",
  2295. digits = c(NA, 2, 0, 4),
  2296. booktabs = TRUE,
  2297. escape = FALSE)%>%
  2298. kableExtra::kable_styling(latex_options = "HOLD_position") %>% row_spec(4,hline_after=TRUE)
  2299. ```
  2301. The experiments were conducted in a windowless room, in which the only source of light was the modified VR headset. The participants sat as close as possible to the VR headset, leaving enough space for the eye-tracker to record the eyes. Each experiment was carried out in a single session of around 45-60 minutes, divided into three blocks of 15-17 minutes each. In each block, there were a total of 60 trials lasting 15 seconds each (12s of stimulus presentation, followed by 3s of interstimulus interval). The participants were given no task other than look at the black fixation dot while trying to minimise their blinking during the presentation period. For all experiments other than the light flux condition, participants adapted to the unmodulated background luminance for two minutes before stimulation began.
  2303. Before the start of the L-M experiment, participants completed a luminance nulling perceptual calibration procedure in L-M cone space on an Iiyama $\textrm{VisionMaster}^{\textrm{TM}}$ Pro 510 display (800 x 600 pixels, 60 Hz refresh rate). During the task, participants were presented with a disc flickering within the L-M cone space (between magenta and cyan). Using a trackball, participants adjusted the angle in cone space to find their subjective isoluminant point, which resulted in changing the flickering intensity of the stimulus until the amplitude of the flicker appeared to be minimised. The result was used to modify the requested contrasts during stimulus preparation so as to account for individual differences affecting perceived illuminance, principally the L:M cone ratio [@Carroll2002; @Hofer2005].
  2305. ## Data analysis
  2307. The pupillometry data were analysed using the same method we used in our previous study [@Segala2023]. The data were converted from mp4 videos to a csv text file using the Pupil Player software [@Kassner2014], which estimated pupil diameter for each eye on each frame using a 3D model of the eyeball. The individual data were then loaded into R for analysis, where a ten-second waveform for each trial in each eye was extracted (excluding the first two seconds after stimulus onset). We interpolated across any dropped or missing frames to ensure regular and continuous sampling over time. The Fourier transform was calculated for each waveform, and all repetitions of each condition were pooled across eye and then averaged. Finally, data were averaged across all participants to obtain the group results. We used coherent averaging and at each stage we excluded data points with a Mahalanobis distance exceeding $D$ = `r SDthresh` from the complex-valued mean [@Baker2021]. For monocular stimulation, we confirmed that the consensual response was equivalent to the response in the stimulated eye.
  2309. For all experiments, we used a bootstrapping procedure with `r nbootstraps` iterations to estimate standard errors across participants. All analysis and figure construction was conducted using a single R-script, available online, making this study fully computationally reproducible: https://osf.io/gdvt4/.
  2311. ## Computational model and parameter estimation
  2313. To quantitatively summarise our data, we used the same model described in our previous study [@Segala2023]. The model has the same general form as the first stage of the contrast gain control model proposed by Meese and colleagues [@Meese2006], but omits the second stage. For the previous model that we used [@Segala2023], the exponent of the numerator and denominator had fixed values of 2 and (implicitly) 1. Here, we allow these parameters (called _p_ and _q_) to be free, to allow different shapes of contrast response function, e.g. accelerating or saturating.
  2315. In words, the model works as follows. Each eye produces an excitatory response that grows with the contrast in that eye. This response is then divided by a term that includes the contrast in the other eye, so that stimulating one eye reduces the response to the other. The two suppressed monocular signals are summed to give the binocular response. The weight $\omega$ sets the strength of this interocular suppression: values below 1 give weak suppression, and therefore stronger binocular facilitation, whereas values near or above 1 give strong suppression, which can produce ocularity invariance or, when suppression is strong enough, a binocular response that is smaller than the monocular one. The responses of the left eye and right eye channels are as follows:
  2317. \begin{equation}
  2318. \label{eq:respL}
  2319. resp_L = \frac{L^p}{Z + L^q + \omega R^q},
  2320. \end{equation}
  2322. \begin{equation}
  2323. \label{eq:respR}
  2324. resp_R = \frac{R^p}{Z + R^q + \omega L^q},
  2325. \end{equation}
  2327. \noindent where _L_ and _R_ are the contrast signals from the left and right eyes, _p_ and _q_ are exponents, _Z_ is a saturation constant that shifts the contrast-response function laterally, and $\omega$ is the weight of suppression from the other eye.
  2329. The responses from the two eyes are then summed binocularly:
  2331. \begin{equation}
  2332. \label{eq:respB}
  2333. resp_B = R_{max}(resp_L + resp_R) + k,
  2334. \end{equation}
  2336. \noindent where _k_ is a noise parameter, and $R_{max}$ scales the overall response amplitude.
  2338. The models were fit using a hierarchical Bayesian framework implemented in Stan [@Carpenter2017]. The data for each photoreceptor type and response frequency was fit separately, for a total of 8 model fits. We did not fit the two eyes separately. Target-eye assignment was counterbalanced across participants, and responses were pooled across eye before fitting, after confirming that the consensual response matched the response in the stimulated eye (see Data analysis). The prior for the $\omega$ parameter was Gaussian, with a mean of 1 and standard deviation of 0.5. Priors for the other free parameters were also Gaussian, with mean values based on previous work [@Segala2023]. All free parameters were constrained to be positive (and the saturation constant _Z_ to be at least 1), so these Gaussian priors were truncated at those bounds. The model was hierarchical and was fitted to the individual-participant amplitudes rather than to the group averages: each participant had their own set of parameters, drawn from group-level distributions whose means we report and compare across conditions. We ran 16 chains of 64,500 iterations each, discarded the first 2,000 iterations of each as warm-up, and thinned by a factor of 10, giving 103,200 retained samples per fit (over `r numSavedSteps*thinSteps` samples before thinning). The group-level parameters converged well in every fit (split-\(\hat{R} < 1.02\) and effective sample size \(> 3{,}000\)). In the melanopsin first-harmonic fit a few participant-level exponent and suppression parameters reached \(\hat{R} = 1.06\); the group-level suppression weight, which is our parameter of interest, had \(\hat{R} = 1.00\). As a summary of fit, we report the root-mean-square error between the group-mean data and the model prediction for each fit in Table \@ref(tab:diagnosticstable).
  2340. ```{r designtable, echo=FALSE}
  2342. design <- data.frame(
  2343. Condition = c("Monocular (0.5 Hz)", "Binocular (0.5 Hz)", "Dichoptic (0.5 Hz)",
  2344. "Monocular (0.4 Hz)", "Binocular cross", "Dichoptic cross"),
  2345. TargetEye = c("0.5 Hz, variable", "0.5 Hz, variable", "0.5 Hz, variable",
  2346. "0.4 Hz, variable", "0.5 Hz, variable", "0.5 Hz, variable"),
  2347. OtherEye = c("Steady background", "0.5 Hz, variable", "0.5 Hz, fixed 48%",
  2348. "Steady background", "0.4 Hz, variable", "0.4 Hz, fixed 48%"),
  2349. AnalysedAt = c("0.5 Hz", "0.5 Hz", "0.5 Hz", "0.4 Hz", "0.5 & 0.4 Hz", "0.5 & 0.4 Hz"),
  2350. Purpose = c("Response to one eye", "Binocular facilitation", "Interocular suppression",
  2351. "Frequency control", "Cross-frequency interaction", "Cross-frequency suppression")
  2352. )
  2354. kable(design,
  2355. caption = "The six ocular conditions used in each of the four experiments. The target eye was modulated at the target contrast (one of five levels); the response was measured at the frequency or frequencies indicated. Target-eye assignment was counterbalanced across the left and right eyes.",
  2356. col.names = c("Condition", "Target eye", "Other eye", "Analysed at", "Purpose"),
  2357. align = "lllll",
  2358. booktabs = TRUE,
  2359. escape = TRUE) %>%
  2360. kableExtra::kable_styling(latex_options = "HOLD_position") %>% row_spec(4,hline_after=TRUE)
  2361. ```
  2364. # Results
  2366. Figure \@ref(fig:spectraplots)m-p shows averaged waveforms of pupil diameter in response to binocular stimulation. For light flux stimuli (Figure \@ref(fig:spectraplots)m) a strong modulation is apparent at the stimulation frequency (0.5Hz), which is also clear in the Fourier amplitude spectrum (Figure \@ref(fig:spectraplots)q). For binocular stimulation of the L-M pathway, S-cone pathway and melanopsin pathway, pupil modulations were less than 10\% of the amplitude of the light flux modulation (note the change in y-axis scale for the Fourier spectra), but still apparent at 0.5Hz in the Fourier spectra (Figure \@ref(fig:spectraplots)n-p). We also observed a second harmonic response (at 1Hz) for all conditions, which was weaker than the first harmonic for light flux and melanopsin stimulation, but stronger for L-M and S-cone stimulation. The second harmonic is also apparent in the pupil waveforms shown in Figure \@ref(fig:spectraplots)n-p. Our main analysis therefore focuses on the amplitude of the pupil modulations at both the first and second harmonic frequencies across different stimulus conditions.
  2368. A response at twice the stimulus frequency tells us that the pupil does not simply follow the sinusoidal modulation. A purely linear system driven by a sinusoid responds only at the stimulus frequency, so energy at the second harmonic must reflect a nonlinear transformation of the input. Such harmonic distortion can arise from many sources including rectification or asymmetry between constriction and redilation. We therefore treat the second harmonic as a complementary measure of nonlinear temporal processing in the pupil response, rather than as a separate oscillatory mechanism. This analysis was exploratory, as we had no specific prediction about the size of the second harmonic in each condition.
  2370. Figure \@ref(fig:FirstHarmonicPlots)a-d shows contrast response functions across stimulation conditions for responses at the first harmonic of the main stimulation frequency (0.5Hz). In each plot, the response to monocular stimulation is given by the red circles and typically increases monotonically as a function of stimulus (temporal) contrast. Relative to monocular stimulation, binocular stimulation led to higher response amplitudes, indicating a binocular facilitation effect, for the light flux and S-cone conditions (Figure \@ref(fig:FirstHarmonicPlots)a,c), and to some extent for the melanopsin condition (Figure \@ref(fig:FirstHarmonicPlots)d). But the L-M cone condition (Figure \@ref(fig:FirstHarmonicPlots)b) produced a binocular suppression effect, where the response to binocular stimulation was weaker than the response to monocular stimulation (blue squares below red circles). These results indicate that the magnitude of binocular facilitation differs across photoreceptor pathway, suggesting heterogeneity in the underlying neural computation.
  2372. ```{r FirstHarmonicPlots, fig.cap="Contrast response functions for pupil modualations in response to flicker at 0.5Hz (panels a-d), and 0.4Hz (panels f-i), and posterior parameter distributions for the weight of interocular suppression (panel e). Within each panel (except panel e), data points are the coherently averaged amplitudes for each condition, and error bars indicate bootstrapped 95\\% confidence intervals. Curves show model fits using the maximum a posteriori (MAP) parameter values (see Table 3). In panel (e), vertical lines show the MAP estimates, and the thin black curve indicates the prior (note the logarithmic x-axis).", fig.align="center", out.width="40%", echo=FALSE}
  2374. knitr::include_graphics(paste0(figdir,'FirstHarmonic.pdf'))
  2376. ```
  2378. In contemporary models of binocular signal combination, the amount of binocular facilitation is determined by the magnitude of interocular suppression, with strong suppression reducing facilitation [@Meese2006]. We can estimate the strength of interocular suppression by measuring how much monocular responses are reduced when a dichoptic 'mask' is shown to the other eye. In our paradigm, the two components flickered at different frequencies (0.5 and 0.4Hz) so that their responses remained distinct in the Fourier spectrum [e.g. @Busse2009] (all conditions are summarised in Table \@ref(tab:diagnosticstable)). The yellow diamond symbols in Figure \@ref(fig:FirstHarmonicPlots)a-d show the target responses in this condition, and in most cases were weaker than the monocular responses (red circles). The strongest dichoptic masking is found in the L-M condition, where we also observed the binocular suppression effect. Suppression can also be estimated from the responses at 0.4Hz (Figure \@ref(fig:FirstHarmonicPlots)f-i). The reduced amplitude in the binocular cross condition (where the two eyes received different temporal frequencies; purple squares) relative to the 0.4Hz monocular condition (grey circles), and the progressive decline in amplitude of the dichoptic cross response (yellow diamonds) also differ across photoreceptor conditions, showing similar differences to those observed at 0.5Hz.
  2380. To estimate the extent of interocular suppression for each photoreceptor pathway, we fitted each data set using a Bayesian hierarchical implementation of a simple binocular combination model [@Meese2006; @Segala2023]. Our primary objective was to compare posterior distributions of the weight of interocular suppression, which are shown in Figure \@ref(fig:FirstHarmonicPlots)e. Consistent with our earlier observations, the strongest suppressive weight corresponds to the L-M and Melanopsin conditions, and the weakest suppression corresponds to the light flux and S-cone conditions, with virtually no overlap between the posterior distributions for weak and strong suppression. The model captured the main features of the contrast response functions (Figure \@ref(fig:FirstHarmonicPlots)a-d,f-i), fitting the light flux data most closely and it fitted the noisier S-cone data least well; it did not capture every data point, and some scatter remains in the weaker responses. Fit statistics and convergence diagnostics for all eight fits are given in TableTable \@ref(tab:diagnosticstable), and a posterior predictive check (Figure \@ref(fig:FS3)) plots the model's 95\% predictive intervals against the data for every condition: most points fall within these intervals, with the widest intervals and the largest deviations in the noisier, low-amplitude conditions. The fitted model parameters are given in Table \@ref(tab:paramtable).
  2382. ```{r paramtable, echo=FALSE}
  2384. condnames <- c('Light flux 1F','L-M cone 1F','S cone 1F','Melanopsin 1F','Light flux 2F','L-M cone 2F','S cone 2F','Melanopsin 2F')
  2385. zlist <- c(pupilZLum, pupilZLM, pupilZS, pupilZMel, pupilZLum2, pupilZLM2, pupilZS2, pupilZMel2)
  2386. klist <- c(pupilKLum, pupilKLM, pupilKS, pupilKMel, pupilKLum2, pupilKLM2, pupilKS2, pupilKMel2)
  2387. wlist <- c(pupilWLum, pupilWLM, pupilWS, pupilWMel, pupilWLum2, pupilWLM2, pupilWS2, pupilWMel2)
  2388. plist <- c(pupilPLum, pupilPLM, pupilPS, pupilPMel, pupilPLum2, pupilPLM2, pupilPS2, pupilPMel2)
  2389. qlist <- c(pupilQLum, pupilQLM, pupilQS, pupilQMel, pupilQLum2, pupilQLM2, pupilQS2, pupilQMel2)
  2390. Rlist <- c(pupilRmaxLum, pupilRmaxLM, pupilRmaxS, pupilRmaxMel, pupilRmaxLum2, pupilRmaxLM2, pupilRmaxS2, pupilRmaxMel2)
  2391. tabledata <- data.frame(condnames,round(zlist,digits=2),round(klist,digits=5),round(wlist,digits=2),round(plist,digits=2),round(qlist,digits=2),round(Rlist,digits=4))
  2393. colnames(tabledata) <- c('Experiment','Z','k','w','p','q','Rmax')
  2394. kable(tabledata, align='lcccccc',booktabs=TRUE,linesep='',caption = 'Summary of maximum a posteriori (MAP) parameter estimates for each data set.') %>%
  2395. kableExtra::kable_styling(latex_options = "HOLD_position") %>% row_spec(4,hline_after=TRUE)
  2397. ```
  2399. At the second harmonic frequencies, the levels of suppression were more uniform across different photoreceptor pathways (see Figure \@ref(fig:SecondHarmonicPlots)). In general, suppression estimates were near or above 1 (see lower rows of Table \@ref(tab:paramtable)), with substantial overlap between the posterior distributions (Figure \@ref(fig:SecondHarmonicPlots)e). The contrast response functions also looked more uniform, and generally involved less binocular facilitation and more interocular suppression than were seen at the first harmonics. The L-M condition now featured the weakest suppressive weight, and the strongest binocular facilitation, which was the opposite pattern seen at the first harmonic.
  2401. Because response energy can shift between the first and second harmonics depending on the exact shape of the transducer, either harmonic on its own could misrepresent the total periodic response. As a check, we combined the two harmonics into a single measure of periodic response magnitude: the root-sum-square of the first- and second-harmonic amplitudes (Figure \@ref(fig:FS2)). The combined measure preserved the main qualitative pattern. Binocular stimulation exceeded monocular stimulation for the light flux and S-cone conditions (binocular facilitation), whereas the L-M and melanopsin conditions showed little facilitation and sat close to ocularity invariance. Combining the harmonics moderated the most extreme ratios, in particular the large S-cone facilitation and the L-M binocular suppression seen at the first harmonic alone, but did not fundamentally alter our conclusions.
  2403. ```{r SecondHarmonicPlots, fig.cap="Responses at the second harmonic frequencies (1Hz and 0.8Hz), in the same format as Figure 2.", fig.align="center", out.width="40%", echo=FALSE}
  2405. knitr::include_graphics(paste0(figdir,'SecondHarmonic.pdf'))
  2407. ```
  2409. Finally, we inspected the response phase of our four stimulation conditions, given previous reports that these differ across pathways [@Spitschan2014], and may be in antiphase for melanopsin and S-cone signals. Figure \@ref(fig:phaseplots) shows the phase angles for the first (a) and second (b) harmonic frequencies for binocular stimulation (monocular stimulation produced very similar results). At the first harmonic, melanopsin and S-cone signals differed in phase by more than 90 degrees, though they were not fully in antiphase. The light flux and L-M responses were approximately in antiphase to each other, however this is likely due to our choice to modulate L+M- in the first half-cycle of the sine wave (corresponding to a luminance increase in the light flux condition), and L-M+ in the second half-cycle (corresponding to a luminance decrease in the light flux condition). Had we reversed this phase arrangement, we would likely have seen a close phase correspondence between the L-M and light flux conditions (another way to think about this is that L-cone decreases and M-cone increases are processed like luminance increases). Phase differences between the light flux, S-cone and melanopsin conditions are likely to reflect, at least in part, different lags in phototransduction at the earliest stage (i.e. in the retina), although downstream temporal filtering within the pupillary pathway could also contribute. At the second harmonic frequency (Figure \@ref(fig:phaseplots)b), the light flux condition was again out of phase with the other three, and was approximately in quadrature phase with the L-M condition, and in antiphase with both the S-cone and melanopsin conditions (which were in phase with each other). We note that the marked phase differences between conditions make the possibility that our results are dominated by rod activity relatively unlikely (e.g. in the L-M and melanopsin conditions, which have very different phase profiles).
  2411. ```{r phaseplots, fig.cap="Pupil phase plots at the first and second harmonic frequencies for the light flux, melanopsin, L-M pathway and the S-cone pathway conditions. Panel (a) shows the pupil response at the first harmonic frequency during binocular stimulation. Panel (b) shows the pupil response at the second harmonic frequency during binocular stimulation. The five data points for each pathway correspond to the five contrast levels displayed. In panel (a) the light flux amplitudes have been scaled down by a factor of 10 to enable comparison with the other conditions.", fig.align="center", echo=FALSE}
  2413. knitr::include_graphics(paste0(figdir,'phaseplots.pdf'))
  2415. ```
  2417. # Discussion
  2419. We used binocular pupillometry and silent substitution to measure monocular and binocular responses of the pupils to flickering stimuli when stimulating specific photoreceptor pathways. In all four experiments, we were able to record contrast response functions at both the first and the second harmonic frequencies. All experiments showed that binocular combination in the autonomic nervous system happens in a non-linear manner, with evidence of different magnitudes of interocular suppression depending on the photoreceptor pathway. This pattern of results was confirmed by a computational model, which allowed us to compare the weight of interocular suppression for each pathway. We found that at the first harmonic frequency the L-M and melanopsin pathways involved strong suppression, whereas the light flux and S-cone pathways involved weaker suppression. Suppression was strong in all four pathways at the second harmonic frequency. Finally, the phase of the pupil response revealed different lag times for the different pathways at the first and second harmonic frequencies.
  2421. This is the first study to investigate binocular interactions in the melanopsin pathway directly. A previous analysis [@Spitschan2019] indicated that there may be a substantial binocular facilitation effect in the circadian pathway, as indexed by melatonin suppression (melatonin is a hormone released by the pineal gland; its production is suppressed by exposure to bright light, particularly when the melanopsin-containing ipRGCs are stimulated). In brief, monocular stimulation of the ipRGCs requires up to ten times the signal strength to produce an equivalent effect to binocular stimulation. This superadditive effect implies an absence of interocular suppression, and perhaps the presence of either a highly compressive nonlinearity, or an AND-style neural operation. Our findings here are different - we do not see substantial binocular facilitation effects in response to melanopsin modulation, and our data indicate that interocular suppression is strong. By measuring the full contrast response function, we can also rule out compressive nonlinearities, because the functions accelerate at both the first and second harmonic frequencies. The explanation for this difference could be due to different anatomical pathways. Binocular combination in the circadian system likely takes place in the suprachiasmatic nucleus, whereas in the pupil constriction circuit the Edinger-Westphal nucleus is the most likely site of binocular integration [@Mathot2018]. Presumably these anatomical differences, and the practical constraints of the two systems, lead to differences in response.
  2423. Our recent psychophysical work [@Baker2024] has looked at binocular interactions in the L-M and S-cone pathways, and compared these to the light flux pathway. Using a contrast discrimination paradigm with spatial modulations of light flux and colour, we found equally strong interocular suppression in all three pathways [@Baker2024]. This is rather different from our pupillometry results here, which show weaker suppression in the light flux and S-cone pathways than in the L-M pathway (see Figure \@ref(fig:FirstHarmonicPlots)). But the current experiments involve temporal modulations, which are quite distinct from the spatial modulations used in our previous work [@Baker2024]. Our other recent work [@Segala2023] has shown that temporal luminance modulations involve much weaker interocular suppression in the cortical response than do spatial luminance modulations [@Baker2017]. These different normalisation processes might reflect different priorities for spatial and temporal vision. Spatial vision aims to fuse images to provide binocular single vision, and benefits from 'ocularity invariance' [@Baker2007], in which visual appearance is constant when viewed with one eye or two. Temporal vision is critical for motion perception, which can involve alternative binocular computations, such as calculating velocity differences between the eyes [@Kaestner2019]. Of course the present measurements were of pupil size, which may be subject to different anatomical and functional constraints from those in the cortex. Future research should aim to extend our current findings to perception using psychophysical approaches.
  2425. ## Physiological interpretation of the pathway and harmonic differences
  2426. Our four experiments modulate different photoreceptor classes, and the differences between them are open to physiological interpretation. We are cautious about attributing responses to discrete physiological mechamisms because silent substitution isolates photoreceptor contrast rather than an anatomically separate neural pathway: the targeted photoreceptor signals converge with others downstream, and we did not explicitly silence the rods although their contribution at the mean luminance levels used is expected to be small.
  2428. The light flux stimulus modulates several photoreceptor classes together and produced by far the largest pupil response. Its weak first-harmonic suppression is consistent with strong binocular facilitation, but because several receptors are co-modulated this condition cannot tell us which retinal input dominates. The isoluminant L-M stimulus targets a cone-opponent signal rather than luminance, and showed the strongest first-harmonic suppression, with a binocular response weaker than the monocular one. This may indicate stronger binocular normalisation for cone-opponent input, though our data cannot localise whether this interaction arises in retinal, pretectal or other subcortical circuitry.
  2430. The S-cone stimulus produced a relatively large second harmonic. This could reflect nonlinearity or asymmetry in how S-cone increments and decrements are converted into a pupil response, rather than a separate oscillatory mechanism. For example, it is well known that there are asymmetries in the number of S-cone on- and off- bipolar cells in the retina [@Klug2003, @LeeGrunert2007]. The melanopsin-directed responses were delayed in phase relative to the cone-directed responses, which is broadly consistent with the slower phototransduction of the ipRGCs. But the ipRGCs also receive cone and rod input, and our stimulus isolates melanopsin contrast rather than an anatomically independent pathway, so the delay need not arise from melanopsin phototransduction alone.
  2432. The meaning of the differences we observe between the 0.4 and 0.5 Hz responses is unclear. The pupil is a temporally filtered and nonlinear system, so a modest change in frequency can alter response phase, waveform asymmetry and harmonic content. Some of the difference may also reflect noisier estimates or imperfect model fit. We therefore do not attribute the 0.4 versus 0.5 Hz differences to distinct physiological mechanisms. A limitation that applies across all of these comparisons is that the four experiments used separate groups of 24 participants, so differences between the photoreceptor-directed conditions are partly between-group differences.
  2434. # Conclusions
  2436. We have shown that binocular combination of temporal flickering light in the autonomic nervous system depends on which photoreceptors are stimulated. We were able to elicit pupil responses by stimulating the periphery of the retina, and to record contrast response functions for all four photoreceptor-directed conditions. While all conditions showed non-linear combination, they varied in how the signals are combined, particularly in the weight of interocular suppression. This was strong ($\omega \ge 1$) for L-M and melanopsin signals at the first harmonic, and all pathways at the second harmonic. Suppression was weaker ($\omega < 1$) in the light flux pathway (consistent with previous work), and also for S-cone directed modulations.
  2438. # References
  2440. <div id='refs'> </div>
  2442. \pagebreak
  2444. # Supplementary material
  2446. \beginsupplement
  2448. ```{r primarylinearity, fig.cap="Linearity of primaries for each device. The calibration measures were summed (i.e., total unweighted irradiance), and the input-output relationship summarised by a 7th order polynomial reverse curve fit. By applying the coefficients of the regression, it is possible to achieve a linear output in the primaries.", fig.align="center", echo=FALSE}
  2450. knitr::include_graphics(paste0(figdir,'gamma_figure.pdf'))
  2452. ```
  2454. ```{r FS2, fig.cap="Combined first- and second-harmonic pupil responses. For each photoreceptor-directed condition, the root-sum-square amplitude of the first (0.5~Hz) and second (1~Hz) harmonic is plotted against target contrast, for monocular, binocular and dichoptic stimulation. Points show the coherently averaged group amplitude and error bars show bootstrapped 95% confidence intervals (10000 iterations). Combining the harmonics preserves the binocular facilitation seen for light flux and S-cone stimulation, and the near-invariance seen for L-M and melanopsin stimulation.", fig.align="center", echo=FALSE}
  2456. knitr::include_graphics(paste0(figdir,'FigureS2.pdf'))
  2458. ```
  2460. ```{r FS3, fig.cap="Posterior predictive check for the eight model fits. Each panel plots the observed group-mean amplitudes (points, with bootstrapped 95% confidence intervals) for the monocular, binocular and dichoptic-cross conditions against target contrast, at the first harmonic (top row) and second harmonic (bottom row) for each photoreceptor-directed condition. Solid lines show the posterior predictive median and shaded regions the 95% posterior predictive interval, computed from the group-level parameters. Most data points fall within the predictive intervals; the intervals are widest, and deviations largest, for the noisier low-amplitude conditions.", fig.align="center", echo=FALSE}
  2462. knitr::include_graphics(paste0(figdir,'FigureS3.pdf'))
  2464. ```

manuscript.Rmd, no license · at the source

Overview

Authors: Federico G Segala1, Aurelio Bruno1,2, Joel T Martin1,3, Anisa Y Morsi1, Alex R Wade1,4, Daniel H Baker1,4
ORCID iDs: Federico G Segala, Aurelio Bruno, Joel T Martin, Anisa Y Morsi, Alex R Wade, Daniel H Baker
  1. Department of Psychology, University of York, York, United Kingdom
  2. School of Psychology and Vision Sciences, University of Leicester, Leicester, United Kingdom
  3. School of Philosophy, Psychology and Language Sciences, University of Edinburgh, Edinburgh, United Kingdom
  4. York Biomedical Research Institute, University of York, York, United Kingdom
Institutions: University of York (United Kingdom); University of Leicester (United Kingdom); University of Edinburgh (United Kingdom)
Journal: PloS one, volume 21, issue 9, article e0356765
Dates: received 18 December 2025; accepted 6 August 2026; published online 1 September 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pone.0356765 · PMID 42679018 · PMCID PMC13533428 · OpenAlex W4399375988
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2024.06.04.597314
Status: code verified
Categories: human (organism)
Methods: Statistics, Physiology & signal measures
MeSH: Autonomic Nervous System*, Vision, Binocular*, Animals, Humans, Photic Stimulation, Pupil, Retinal Ganglion Cells (* major topic)
Topic: Glaucoma and retinal disorders (Ophthalmology, Medicine), according to OpenAlex
Funding: Biotechnology and Biological Sciences Research Council (BB/V007580/1)
Citations: not cited yet (Europe PMC); 52 references in the paper

Abstract

Pupil diameters are regulated by the autonomic nervous system, which combines light signals across the eyes independently of the visual cortex. Distinct classes of retinal photoreceptor are involved in this process, with cones and rods driving the initial constriction and intrinsically photosensitive retinal ganglion cells maintaining diameter over prolonged time periods. We investigated binocular combination by targeting different photoreceptor pathways using a novel binocular multiprimary system to modulate the input spectra via silent substitution. At the first harmonic of the modulation frequency, luminance and S-cone responses showed strong binocular facilitation, and weak interocular suppression. Melanopsin responses were invariant to the number of eyes stimulated. The L-M pathway involved binocular inhibition, whereby responses to binocular stimulation were weaker than for monocular stimulation. The second harmonic involved strong interocular suppression in all pathways, but with some evidence of binocular facilitation. Our results are consistent with a computational model of binocular signal combination (implemented in a Bayesian hierarchical framework), in which the weight of interocular suppression differs across pathways. We also find pathway differences in response phase, consistent with different lag times for phototransduction. These results provide evidence that binocular interactions in the pupillary pathway differ across photoreceptor-directed modulations and across harmonic components of the response.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 16 matches between paragraphs and lines of code.

OSF gdvt4

License: none: the authors keep all their rights
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Languages: Python (9), R (1)
Size: 56 files, 10 scripts
Software Heritage: not checked
Found in: the text, “Data analysis”
Holds: README, license file, environment (docker-compose.yml, docker/Dockerfile), continuous integration, 1 notebook
Not found: CITATION.cff, tests, documentation
Tools: pandas (9 files), Matplotlib (6 files), NumPy (5 files), SciPy (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers (HTTP 200)
  • 26 September 2026: the link answers (HTTP 200)
12 files
At the source: osf.io/gdvt4/

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 10 scripts, each with its path and the digest of its content;
  • 16 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

All relevant data is available on OSF at https://www.doi.org/10.17605/OSF.IO/GDVT4.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 3, 28 September 2026

  • Authors: added Joel T Martin (0000-0002-4475-3835); removed Joel T Martin

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 7 MeSH terms, 1 funder, 51 references.

Cite

This paper

Segala, F. G., Bruno, A., Martin, J. T., Morsi, A. Y., Wade, A. R., & Baker, D. H. (2026). Binocular combination in the autonomic nervous system. PloS one, 21(9), e0356765. https://doi.org/10.1371/journal.pone.0356765

BibTeX

@article{segala2026binocular,
author = {Segala, Federico G and Bruno, Aurelio and Martin, Joel T and Morsi, Anisa Y and Wade, Alex R and Baker, Daniel H},
title = {{Binocular combination in the autonomic nervous system}},
journal = {PloS one},
year = {2026},
month = sep,
volume = {21},
number = {9},
pages = {e0356765},
publisher = {PLOS},
issn = {1932-6203},
doi = {10.1371/journal.pone.0356765},
url = {https://doi.org/10.1371/journal.pone.0356765},
pmid = {42679018},
pmcid = {PMC13533428}
}

RIS

TY - JOUR
AU - Segala, Federico G
AU - Bruno, Aurelio
AU - Martin, Joel T
AU - Morsi, Anisa Y
AU - Wade, Alex R
AU - Baker, Daniel H
TI - Binocular combination in the autonomic nervous system
T2 - PloS one
J2 - PLoS One
PY - 2026
DA - 2026/09/01
VL - 21
IS - 9
SP - e0356765
SN - 1932-6203
PB - PLOS
DO - 10.1371/journal.pone.0356765
UR - https://doi.org/10.1371/journal.pone.0356765
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pone.0356765",
"type": "article-journal",
"title": "Binocular combination in the autonomic nervous system",
"container-title": "PloS one",
"author": [
{
"family": "Segala",
"given": "Federico G"
},
{
"family": "Bruno",
"given": "Aurelio"
},
{
"family": "Martin",
"given": "Joel T"
},
{
"family": "Morsi",
"given": "Anisa Y"
},
{
"family": "Wade",
"given": "Alex R"
},
{
"family": "Baker",
"given": "Daniel H"
}
],
"container-title-short": "PLoS One",
"volume": "21",
"issue": "9",
"page": "e0356765",
"DOI": "10.1371/journal.pone.0356765",
"PMID": "42679018",
"PMCID": "PMC13533428",
"ISSN": "1932-6203",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pone.0356765",
"language": "en",
"issued": {
"date-parts": [
[
2026,
9,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1371/journal.pone.0350140 [code]
Signal combination in flutter vibration perception.
Journal: PloS one
In common: pandas, SciPy, Matplotlib, 1 other tool, 8 references, 2 authors
[2] doi:10.2147/opth.s590186 [code]
Differential Effects of Balanced and Imbalanced Binocular Stimulation on Visual Cortex Responses in Amblyopic Children.
Journal: Clinical ophthalmology (Auckland, N.Z.)
In common: pandas, SciPy, Matplotlib, 1 other tool, 2 references
[3] doi:10.1073/pnas.2538139123
Cognition does not automatically influence perception: Evidence from neural encoding of colors belonging to different categories.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: author Joel T Martin
[4] doi:10.3390/bioengineering13070820 [code]
Benchmarking Multimodal Workload Classification: Effects of Modality, Validation Protocol, and Segmentation Contrast on an Open Graded-Arithmetic Dataset.
Journal: Bioengineering (Basel, Switzerland)
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference
[5] doi:10.7554/elife.108223 [code]
Two time scales of adaptation in human learning rates.
Journal: eLife
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference
[6] doi:10.3758/s13428-026-03122-w [code]
Modeling intraindividual variability in affect (MIVA): Formalized theoretical approach, computational model, and parameter recovery study.
Journal: Behavior research methods
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference
[7] doi:10.1093/braincomms/fcag218 [code]
Magnetoencephalography biomarkers for assessing myelin content and neuronal function in acute optic neuritis.
Journal: Brain communications
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference
[8] doi:10.1126/sciadv.aea3919 [code]
Hierarchical brain dynamics supporting visual perceptual transitions.
Journal: Science advances
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference
[9] doi:10.7554/elife.105953 [code]
Top-down feedback in deep neural networks leads to functional differences during audiovisual integration.
Journal: eLife
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference
[10] doi:10.1038/s41467-026-76098-y [code]
A single computational objective can produce specialization of streams in visual cortex.
Journal: Nature communications
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.