Introduction

Previous work from our lab has shown that saccadic burst generators can have abnormal preferred directions in monkey models of infantile pattern strabismus. Additionally, microstimulation of the superior colliculus and saccadic burst generators can evoke disconjugate movements. This, along with behavioral evidence, suggests a neural origin of the “cross-talk” seen in pattern strabismus.

The intersitial nucleus of cajal (INC) is the location of the vertical neural integrator for eye position. If abnormal signals are present in the burst generators, then these signals should be integrated by the INC. In this report, we assess the direction preference of INC neurons in strabismus and compare with the normal monkey.

library(ggplot2)
library(dplyr)
library(knitr)
library(tidyr)
library(broom)
# library(grid)
library(relaimpo)
library(leaps)
library(stringr)
library(cladoRcpp)
# library(boot)
source('Adamhelperfunctions.R')
select<- dplyr::select

#skip the next two chunks (quickload and mark saccades) if you load this.
 # t<- readRDS('INCsaccadesmarked.RDS')
 t<- readRDS('INCsaccadesmarkedBI2.RDS')

LoadFromFile<- FALSE

if (LoadFromFile){
  o<- readRDS('INCstrab.RDS')
}else{
  o<-NULL
}

#loads all the .csv files in the specified folder, assigns a neuron name to each based 
#on the file name and then binds them all into one data frame. 
#if a reference file is given, then it only loads new data files

t<- loadnewcsv2(path="C:/Users/setup/Desktop/NRTP Vergence/INCstrab/",referencefile=o)

t<- rbind(o,t)

o<-NULL

measureSaccadesINC<-function(t){
  # I set up a check because I want to apply blink removal only for data files 
  #that come from visual eye tracking. There were errors when I tried to run normal 
  #data through it, possibly since it didn't find any periods where the eye postiion 
  #matched the criteria for a blink. I could edit my markSaccades function to handle
  #what to do if there are no saccades found in the entire file, or make a separate 
  #"markblinks" function, but this works well.
  if (t$monkey[1]=='DC'){ 
    t%>%
      mutate(blinks=markSaccades(rep,buffer=80,threshold=70), 
             #we put an 80ms buffer around each detected blink
             #then we force the velocity to be above threshold for the detected blinks
             #This causes all blinks to be marked as saccades, which means they are not
             #included in any fixation analyses
             #When analyzing blinks, remember to omit saccades where the eye postions are extreme
             dsnum=markSaccadesDouble(replace(conj.velocity,blinks>0,150),
                                      threshold1=30,threshold2=15,maxreject=100000),
             sdf=spikedensity(rasters,10),
             sdf10=lag(sdf,10),
             conj.vert=(repV+lepV)/2,
             conj.vert.vel=(revV+levV)/2,
             conj.hor=(rep+lep)/2,
             conj.hor.v=(rev+lev)/2) %>%
      select(-blinks)->
      t
  }else{
    t %>%
      mutate(dsnum=markSaccadesDouble(conj.velocity,threshold1=30,threshold2=15),
              sdf=spikedensity(rasters,10),
              sdf10=lag(sdf,10),
              conj.vert=(repV+lepV)/2,
              conj.vert.vel=(revV+levV)/2,
              conj.hor=(rep+lep)/2,
              conj.hor.v=(rev+lev)/2) ->
      t
  }
  t
}

expandSaccades<- function(tt,buffer=80){
  #this little function adds a new grouping variable
  #so that we can analyze what happens after saccades
  #without losing our ability to just grab the saccades with group_by
  #This way we can group_by(dsnum) to get just the saccades
  #or group_by(dsnum_extended) to get the saccade plus some period after 
  
  maxtime<- max(tt$time)
  tt %>%
    filter(dsnum>0) %>%
    group_by(dsnum) %>%
    summarize(saccade.end=last(time),
              saccade.start=first(time)) %>%
    filter(saccade.end+buffer<maxtime) %>%
    mutate(saccade.end=saccade.end+80)-> 
    tcount
  
tcount$saccade.end<-tcount$saccade.end+buffer

jsac<- function(stimes){
  df<- data.frame(time=stimes[['saccade.start']]:stimes[['saccade.end']],
                  dsnum_extended=stimes[['dsnum']])
  return(df)
}

x<-rbindlist(apply(tcount,1,jsac))

tt<- left_join(tt,x,by='time')

}

#The next block took 2 hours to run
#I needed to increase the memory.limit() otherwise it would error
#I was able to run it by just doubling the limit with  memory.limit(32652)
#Doing this causes the whole computer to freeze up for the two hours while it calculate
#I'm not sure how much time was added by Sophos simultanesouly trying to virus scan
#and the WDBackupEngine doing whatever it does. If I had to run this again I would 
#be sure to disable both of those programs.
#It also took over an hour to simply run saveRDS(t,''), so I think there were some
#issues with memory management in executing this code. 
t %>%
  group_by(neuron) %>%
  do(measureSaccadesINC(.)) %>%
  do(expandSaccades(.,buffer=80))->
  t

# saveRDS(t,'INCsaccadesmarked.RDS')

t %>%
  group_by(neuron)%>%
  mutate(time=row_number()) %>%
  filter(dsnum<0) %>% #fixations only
  group_by(neuron,dsnum) %>%
  summarize(sd.conj.velocity=sd(conj.velocity),
            mean.conj.velocity=mean(conj.velocity),
            spread=max(conj.velocity)-min(conj.velocity),
            qrange=quantile(conj.velocity,0.975)-quantile(conj.velocity,0.025),
            dur=n(),
            starttime=first(time),
            c2eyes=cor(repV,lepV),
            meanFR=sum(rasters)/dur*1000,
            mean.V=mean(conj.vert),
            mean.H=mean(conj.hor),
            mean.R.H=mean(rep),
            mean.L.H=mean(lep),
            mean.R.V=mean(repV),
            mean.L.V=mean(lepV),
            mean.T.V=mean(tvp),
            mean.T.H=mean(thp),
            sd.sdf=sd(sdf10),
            asleep=sd.conj.velocity>7.5 || dur>2000)->
  zp

#In this chunk, we calculate direction preferences for each  neuron based on the measured fixations
conj.measureCell<- function(x,maxEP=25){
  #Remove bad fixations that are too short or are outside of the normal range of eye movements
  # x<- filter(x,dur>200,abs(mean.R.H)<50,abs(mean.L.H)<50)
  x<- filter(x,dur>200,abs(mean.R.H)<maxEP,abs(mean.L.H)<maxEP)
  
  #calculate using mean eye position:
  mH<-lm(meanFR~mean.H,data=x)
  mV<-lm(meanFR~mean.V,data=x)
  mHV<-lm(meanFR~mean.H+mean.R.V,data=x)
  
  H.p=anova(mH)$'Pr(>F)'[1]
  V.p<- anova(mV)$'Pr(>F)'[1]
  H<-as.numeric(coef(mHV)[2])
  V<-as.numeric(coef(mHV)[3]) 
  
  #calculate direction preference based on relative slopes from multiple regression model
  dir.pref.slope<-as.numeric(atan2(V,H))*180/pi

  #calculate direction preference based on relative importance of each factor in multiple regression model  
  b=as.numeric(calc.relimp(mHV)$lmg)
  dir.pref.imp<- atan2(b[2]*sign(V),b[1]*sign(H))*180/pi
  
  any.sig=V.p<0.001 | H.p <0.001
  
  R2<- summary(mHV)$r.squared
  
  H.R2<- summary(mH)$r.squared #conj horizontal fit only
  V.R2<- summary(mV)$r.squared #conj vertical fit only
  
  d<-data.frame(neuron=unique(x$neuron),
                any.sig,
                H.p,V.p,H,V,dir.pref.imp,dir.pref.slope,
                R2,H.R2,V.R2)
                
}

measureCell<-function(x){
  #Remove bad fixations that are too short or are outside of the normal range of eye movements
  x<- filter(x,dur>200,abs(mean.R.H)<50,abs(mean.L.H)<50)
  
  #calculate linear regression models on horizontal, vertical and combined
  #calculate each eye separately 
  mRH<-lm(meanFR~mean.R.H,data=x)
  mLH<-lm(meanFR~mean.L.H,data=x)
  mRV<-lm(meanFR~mean.R.V,data=x)
  mLV<-lm(meanFR~mean.L.V,data=x)
  mRHV<-lm(meanFR~mean.R.H+mean.R.V,data=x)
  mLHV<-lm(meanFR~mean.L.H+mean.L.V,data=x)
  
  #calculate using mean eye position:
  mH<-lm(meanFR~mean.H,data=x)
  mV<-lm(meanFR~mean.V,data=x)
  mHV<-lm(meanFR~mean.H+mean.R.V,data=x)

  
  #Summarize desired information about models
  RH.p=anova(mRH)$'Pr(>F)'[1]
  LH.p=anova(mLH)$'Pr(>F)'[1]
  LV.p<- anova(mLV)$'Pr(>F)'[1]
  RV.p<- anova(mRV)$'Pr(>F)'[1]
  RH<-as.numeric(coef(mRHV)[2])
  RV<-as.numeric(coef(mRHV)[3]) 
  LH<-as.numeric(coef(mLHV)[2])
  LV<-as.numeric(coef(mLHV)[3])  
  
  H.p=anova(mH)$'Pr(>F)'[1]
  V.p<- anova(mV)$'Pr(>F)'[1]
  H<-as.numeric(coef(mHV)[2])
  V<-as.numeric(coef(mHV)[3]) 

  
  
  #two methods to calculate direction preferences:

  #calculate direction preference using slopes of regression lines
  dir.pref.R<-as.numeric(atan2(RV,RH))*180/pi
  dir.pref.L<-as.numeric(atan2(LV,LH))*180/pi
  
  #calculate direction preference using relative importance of horizontal and vertical
  br=as.numeric(calc.relimp(mRHV)$lmg)
  bl=as.numeric(calc.relimp(mLHV)$lmg)
  dir.imp.R<- atan2(br[2]*sign(RV),br[1]*sign(RH))*180/pi
  dir.imp.L<- atan2(bl[2]*sign(LV),bl[1]*sign(LH))*180/pi
  
  any.sig=RV.p<0.001 | LV.p<0.001| LH.p <0.001| RH.p<0.001
  R2R<-summary(mRHV)$r.squared
  R2L<- summary(mLHV)$r.squared
  R2<- summary(mHV)$r.squared
  H.R2<- summary(mH)$r.squared #conj horizontal fit only
  V.R2<- summary(mV)$r.squared #conj vertical fit only
  
  #Calculate RMSE - compare with R-squared?
  # rmseR=sqrt(c(crossprod(mRHV$residuals))/mRHV$df.residual)
  # rmseL=sqrt(c(crossprod(mLHV$residuals))/mLHV$df.residual)
  rmseR=summary(mRHV)$sigma
  rmseL=summary(mLHV)$sigma
  
  #Calculate relationship between firing rate and standard deviation
  msd=lm(sd.sdf~meanFR,data=x)
  sd.int=as.numeric(coef(msd)[1])
  sd.slope=as.numeric(coef(msd)[2])
  
  #This data frame that I am returning is getting kind of big. 
  #I plan to remove the unused measurements once I am past the exploratory phase
  d<-data.frame(neuron=unique(x$neuron),
                RH=RH,RV=RV,LH=LH,LV=LV,dir.pref.R=dir.pref.R,dir.pref.L=dir.pref.L,
                RH.p=RH.p,LH.p=LH.p,RV.p=RV.p,LV.p=LV.p,any.sig=any.sig,
                H.p=H.p,V.p=V.p,H=H,V=V,
                R2R,R2L,dir.imp.R,dir.imp.L,R2=R2,
                rmseR=rmseR,rmseL=rmseL,sd.int,sd.slope)
}

zp %>%
  group_by(neuron) %>%
  do(measureCell(.)) %>%
  separate(neuron,c('monkey','cellnum'),remove=FALSE)->
  zm

zp %>%
  group_by(neuron) %>%
  do(conj.measureCell(.)) %>%
  separate(neuron,c('monkey','cellnum'),remove=FALSE)->
  zmc

Results

Fixation Analysis

A significant challenge in comparing the activity of INC neurons between normal and strabismic animals is the variety of neuronal types present in the normal INC. So-called irregular cells have been identified with lower firing rates and less consistent firing rates. Some of these cells may be burst neurons that become tonically active when the animals are drowsy. Others may have activity related to the position or movement of the head. Fukushima (1990a) described pitch cells with higher coefficients of variation (standard deviation divided by mean of firing rate) and lower mean firing rates in cats. Because our animals did not move their heads.

Note about regression

Even though we are using a simple linear regression analysis on this dataset, we still must make several choices regarding the regression. For each cell, we can compare the firing rate to the Right eye, the Left eye or the cyclopean (mean) eye position. Further, we can regress horizontal and vertical using individual lines (FR ~ H or FR ~ V), or we can use multiple regression to fit the plane that includes both horizontal and vertical eye position (FR ~ H + V).

In the following plot, notice that when we plot firing rate as a function of horizontal eye position for the first cell from the normal animal, the firing rate actually reflects the vertical eye position, which makes the plot reflect the radial spread of target locations. If there is any bias in the presentation of targets, it will appear that the cell has a sensitivity to horizontal eye position. This is especially important for analyzing the data from the strabismic animals as they did not always look equally at targets in all locations. For this reason, considering both horizontal and vertical eye position at the same time is an essential test for the cell’s preferred direction.

Regular Cells

Ncell<- 'Bee-106'
Scell<- 'DC-902'

maxEP<- 25
zp<- filter(zp,
            dur>200,
            abs(mean.R.H)<maxEP,
            abs(mean.L.H)<maxEP,
            abs(mean.R.V)<maxEP,
            abs(mean.L.V)<maxEP,
            meanFR<500)

tab<-zmc %>%
  filter(neuron %in% c(Ncell,Scell)) %>%
  select(neuron,R2,H.R2,V.R2,dir.pref.slope,dir.pref.imp)

multiplot(plotlist=compare2cells(zp,zmc,Ncell,Scell,color1='purple',color2='orange',
                                 c1name='Normal',c2name='Exotrope'),cols=2)

kable(tab,digits=2)

neuron	R2	H.R2	V.R2	dir.pref.slope	dir.pref.imp
Bee-106	0.83	0.01	0.83	-87.90	-89.79
DC-902	0.54	0.04	0.59	118.84	98.42

Above are examples of regular vertical burst tonic cells from a normal animal and an exotrope.

Irregular Cells

Ncell<- 'Bee-118'
Scell<- 'DC-910'

tab<-zmc %>%
  filter(neuron %in% c(Ncell,Scell)) %>%
  select(neuron,R2,H.R2,V.R2,dir.pref.slope,dir.pref.imp)

multiplot(plotlist=compare2cells(zp,zmc,Ncell,Scell,color1='purple',color2='orange',
                                 c1name='Normal',c2name='Exotrope'),cols=2)

kable(tab,digits=2)

neuron	R2	H.R2	V.R2	dir.pref.slope	dir.pref.imp
Bee-118	0.29	0.07	0.24	53.37	75.34
DC-910	0.45	0.11	0.39	-102.57	-98.90

Direction Preferences

Now let’s plot the direction preferences for all the cells using this method:

ggplot(zmc)+
  geom_segment(aes(y=0,yend=1,x=dir.pref.imp,xend=dir.pref.imp,color=monkey))+
  scale_x_continuous(limits=c(-180,180),breaks=c(0,90,180,-90))+
  coord_polar(start=pi/2,direction=-1)+
  theme_minimal()+
  theme(legend.position='none')+
  facet_wrap(~monkey,ncol=1)+
  xlab('Direction preference\nusing relative importance')

But some of these cells do not have a detectable sensitivity to eye position. Using a threshold of p<0.001, we remove 18/83 of these cells:

ggplot(zmc %>% filter(any.sig))+
  geom_segment(aes(y=0,yend=1,x=dir.pref.imp,xend=dir.pref.imp,color=monkey))+
  scale_x_continuous(limits=c(-180,180),breaks=c(0,90,180,-90))+
  coord_polar(start=pi/2,direction=-1)+
  theme_minimal()+
  theme(legend.position='none')+
  facet_wrap(~monkey,ncol=1)+
  xlab('Direction preference\nusing relative importance')

There appears to be a negative correlation between abnormal direction preference and goodness of fit.

ggplot(zmc %>% filter(any.sig))+
  geom_point(aes(90-abs(dir.pref.imp),R2,color=monkey),size=2)+
  theme_minimal()+
  xlab('Direction preference away from vertical')

Let’s scale these lines by the goodness of fit:

ggplot(zmc %>% filter(any.sig))+
  geom_segment(aes(y=0,yend=R2,x=dir.pref.imp,xend=dir.pref.imp,color=monkey))+
  scale_x_continuous(limits=c(-180,180),breaks=c(0,90,180,-90))+
  coord_polar(start=pi/2,direction=-1)+
  theme_minimal()+
  theme(legend.position='none')+
  facet_wrap(~monkey,ncol=1)+
  xlab('Direction preference\nusing relative importance')

Abnormal Direction Preferences

Let’s look at some of the cells with abnormal direction preferences. First I’ll make the same plots as above. Then I will look into the differences in interspike intervals.

W1<-'Bee-118'
W2<-'DC-915'


tab<-zmc %>%
  filter(neuron %in% c(W1,W2)) %>%
  select(neuron,R2,H.R2,V.R2,dir.pref.slope,dir.pref.imp)

multiplot(plotlist=compare2cells(zp,zmc,W1,W2,color1='purple',color2='orange',
                                 c1name='Normal',c2name='Exotrope'),cols=2)

kable(tab,digits=2)

neuron	R2	H.R2	V.R2	dir.pref.slope	dir.pref.imp
Bee-118	0.29	0.07	0.24	53.37	75.34
DC-915	0.17	0.17	0.04	-173.32	-166.94

Evaluating Irregular Cells - CV plots

What is the source of this additional variabilty? In the normal animals, “pitch cells” have been described that respond to movements of the head, and “irregular” cells that seem to not encode spontaneous eye movements but do respond to visual tracking. It’s difficult to say how those cell types might related to the strabismic animal simply because I don’t have recordings of those cell types from a normal monkey using the same setup to be able to compare. I can only look at the figures and tables plotted in published works, but those data are limited.

We can, however, compare the cells from our three monkeys using the criteria described by Fukushima. He discovered that by using the coefficient of variation (CV) and the mean firing rate around the primary position of the eyes. The CV is the standard deviation of the firing rate divided by the mean firing rate. A higher CV means that the cell fires less consistently and is associated with irregular or “pitch” cells in the cat.

plotCells<- c('Bee-106','DC-902','Bee-118','DC-915')
toPlot<-list()
for (i in seq_along(plotCells)){
  zp %>%
    filter(neuron==plotCells[i],meanFR>1) %>%
    ggplot()+
    geom_point(aes(mean.V,sd.sdf/meanFR))+
    ylab('CV=SD/FR')+
    ggtitle(plotCells[i])->
    toPlot[[i]]
}

multiplot(plotlist=toPlot,cols=2)

zp %>% #just look at fixations within 5 deg of zero
  filter(neuron %in% plotCells,meanFR>1,abs(mean.V)<5) ->
  forCVplot

forCVplot %>%
  group_by(neuron) %>%
  summarize(mean.CV=mean(sd.sdf/meanFR,na.rm=T),
            mean.FR=mean(meanFR,na.rm=T))->
  cvmeans

forCVplot %>%
  ggplot()+
  geom_point(aes(mean.V,sd.sdf/meanFR))+
  facet_wrap(~neuron,ncol=2)+
  geom_hline(aes(yintercept = mean.CV),data=cvmeans)+
  ylab('CV=SD/FR')+
  theme_minimal()

kable(cvmeans,digits = 2)

neuron	mean.CV	mean.FR
Bee-106	0.42	27.63
Bee-118	0.41	36.17
DC-902	0.21	49.94
DC-915	0.39	43.30

This figure from Fukushima et al 1991. These data are from cats.

FukushimaCatCV

zp %>%
  filter(meanFR>1,abs(mean.V)<5) %>%
  separate(neuron,c('monkey','cellnum'),remove=FALSE) %>%
  group_by(neuron,monkey) %>%
  summarize(mean.CV=mean(sd.sdf/meanFR,na.rm=T),
            mean.FR=mean(meanFR,na.rm=T))->
  cvmeans

zmc %>%
  select(neuron,R2) %>%
  left_join(cvmeans,.,by='neuron')->
  cvmeans

ggplot(cvmeans)+
  geom_point(aes(mean.CV,mean.FR,color=monkey,size=R2),alpha=0.5)+
  xlim(0,2)+
  ylim(0,150)+
  scale_x_log10()

# cvmeans %>%
#   filter(mean.CV>0) %>%
#   arrange(mean.CV) %>%
#   select(-monkey) %>%
#   kable(digits=2)

To examine what is happening here, let’s look at the two cells from the normal animal that have the highest and lowest mean CV.

multiplot(plotlist=compare2cells(zp,zmc,'Bee-123','Bee-124',color1='purple',color2='orange'),cols=2)

As you can see from this plot, the explanation for the extremely low firing rate and high CV for the first cell is that the cell does not begin to encode vertical eye position until the eyes are above the central position. Why weren’t any cells of this type described in previous papers?

We will have to take x-intercept into account for this analysis. Do the previous studies define the “primary position in the orbits” differently depending on the firing rate of the cells?

low CV cells in strabismus

multiplot(plotlist=compare2cells(zp,zmc,'Kopachuck-907','Kopachuck-943',color1='purple',color2='orange'),cols=2)

CV outside of primary position

Next, I want to look at the CV plot for something other than the “primary” position of the eyes. I’m thinking to either show -10, 0, +10, or do a calculation based on the x intercept. If slope > 0 and x_intercept > 0, then use +10. If slope < 0 and x_intercept < 0 then use -10. Otherwise use zero. Use zero if \(sign(slope)*sign(x_intercept) < 0\), otherwise use \(10*sign(slope)\)

zp %>%
  filter(meanFR>1,abs(mean.V)<5) %>%
  separate(neuron,c('monkey','cellnum'),remove=FALSE) %>%
  group_by(neuron,monkey) %>%
  summarize(mean.CV=mean(sd.sdf/meanFR,na.rm=T),
            mean.FR=mean(meanFR,na.rm=T))->
  cvmeans

zmc %>%
  select(neuron,R2) %>%
  left_join(cvmeans,.,by='neuron')->
  cvmeans

zp %>%
  filter(meanFR>1,mean.V > 8, mean.V< 12) %>%
  separate(neuron,c('monkey','cellnum'),remove=FALSE) %>%
  group_by(neuron,monkey) %>%
  summarize(mean.CVpos=mean(sd.sdf/meanFR,na.rm=T),
            mean.FRpos=mean(meanFR,na.rm=T))->
  cvmeanspositive

zp %>%
  filter(meanFR>1,mean.V> -12,mean.V< -8) %>%
  separate(neuron,c('monkey','cellnum'),remove=FALSE) %>%
  group_by(neuron,monkey) %>%
  summarize(mean.CVneg=mean(sd.sdf/meanFR,na.rm=T),
            mean.FRneg=mean(meanFR,na.rm=T))->
  cvmeansnegative

cvmeans<- left_join(cvmeans,cvmeanspositive,by=c('neuron','monkey'))
cvmeans<- left_join(cvmeans,cvmeansnegative,by=c('neuron','monkey'))

cvmeans %>%
  mutate(mean.CVmin=min(mean.CV,mean.CVpos,mean.CVneg),
         mean.FRmax=max(mean.FR,mean.FRpos,mean.FRneg),
         upward= mean.FRpos>mean.FRneg)->
  cvmeans

p1<- qplot(mean.CV,mean.FR,data=cvmeans,color=monkey)+ggtitle('around zero')+
  scale_x_log10()+geom_vline(xintercept = 0.5)+theme(legend.position='bottom')

p2<- qplot(mean.CVpos,mean.FRpos,data=cvmeans,color=monkey)+ggtitle('up')+
  scale_x_log10()+geom_vline(xintercept = 0.5)+theme(legend.position='bottom') 

p3<- qplot(mean.CVneg,mean.FRneg,data=cvmeans,color=monkey)+ggtitle('down')+
  scale_x_log10()+geom_vline(xintercept = 0.5)+theme(legend.position='bottom')

p4<- qplot(mean.CVmin,mean.FRmax,data=cvmeans,color=monkey)+ggtitle('min')+
  scale_x_log10()+geom_vline(xintercept = 0.5)+theme(legend.position='bottom')

multiplot(plotlist=list(p1,p2,p3,p4),cols = 2)

qplot(mean.CVmin,mean.FRmax,data=cvmeans,color=monkey,size=R2)+ggtitle('min')+
  scale_x_log10()+geom_vline(xintercept = 0.5)+theme(legend.position='bottom')

So Doing this adjustment puts almost all of the cells from the normal animal below the threshold identified by Fukushima in cats for the “pitch cells,” but there is one cell that is still above it. Let’s plot it now and compare it to one of the still-irregular cells from the Esotrope.

multiplot(plotlist=compare2cells(zp,zmc,'Bee-101','Kopachuck-943',color1='purple',color2='orange'),cols=2)

Saccade Analysis

The previous analysis was based on the mean firing rates of the cells while the eyes were held in various positions in the orbits. In the next section, we will look at the dynamic activity of the cells during saccades. The deeper analyses in this section will ncessarily be limited to the subset of cells who respond robustly with a burst of spikes during saccades. To identify these cells, let’s first plot the peak firing rate as a function of peak gaze velocity, split into horizontal and vertical components.

Peak Firing Rate Analysis

plotsaccades<-function(bvelplot,cellname,plotColor='blue'){
  multiplot(
    bvelplot%>%
      filter(neuron==cellname)%>%
      ggplot()+
      geom_point(aes(peak.V.Velocity,avgFR),alpha=0.8,color=plotColor)+
      stat_smooth(aes(peak.V.Velocity,avgFR,group=neuron),color='black',method='lm',se = FALSE)+
      ylim(0,NA)+
      xlab('Peak Vertical Eye Velocity (deg/s)')+
      ylab('Firing Rate during Saccade\n(spk/s)')+
      ggtitle(cellname)
    ,
    bvelplot%>%
      filter(neuron==cellname)%>%
      ggplot()+
      geom_point(aes(peak.H.Velocity,avgFR),alpha=0.8,color=plotColor)+
      stat_smooth(aes(peak.H.Velocity,avgFR,group=neuron),color='black',method='lm',se = FALSE)+
      ylim(0,NA)+
      xlab('Peak Horizontal Eye Velocity (deg/s)')+
      ylab('Firing Rate during Saccade/n(spk/s)')
  )
}
# t %>%
#   group_by(neuron)%>%
#   mutate(time=row_number(),
#          lagspikes=lag(rasters,10)) %>%
#   filter(dsnum>0) %>% #saccades only
#   group_by(neuron,dsnum) %>%
#   summarize(peak.conj.velocity=maxabs(conj.velocity),
#             dur=n(),
#             starttime=first(time),
#             peakFR=max(sdf10),
#             nspikes=sum(lagspikes),
#             #for the eye tracker, whenever there is a blink it puts in values > 60. Reject those.
#             maxep=max(abs(rep),abs(lep),abs(repV),abs(lepV)))->
#   zp


# t %>%
#   group_by(neuron)%>%
#   mutate(time=row_number()) %>%
#   filter(dsnum>0) %>% #saccades
#   group_by(neuron,dsnum) %>%
#   summarize(sd.conj.velocity=sd(conj.velocity),
#             mean.conj.velocity=mean(conj.velocity),
#             sd.verg.velocity=sd(verg.velocity),
#             spread=max(conj.velocity)-min(conj.velocity),
#             qrange=quantile(conj.velocity,0.975)-quantile(conj.velocity,0.025),
#             dur=n(),
#             R.H.Amp=last(rep)-first(rep),
#             R.V.Amp=last(repV)-first(repV),
#             L.H.Amp=last(lep)-first(lep),
#             L.V.Amp=last(lepV)-first(lepV),
#             disjunctiveH=sign(R.H.Amp*L.H.Amp)<0,
#             disjunctiveV=sign(R.V.Amp*L.V.Amp)<0,
#             conj.H.Amp=(R.H.Amp+L.H.Amp)/2,
#             conj.V.Amp=(R.V.Amp+L.V.Amp)/2,
#             peak.conj.velocity=maxabs(conj.velocity),
#             peak.H.Velocity=maxabs((rev+lev)/2),
#             peak.V.Velocity=maxabs((revV+levV)/2),
#             peak.RH.Velocity=maxabs(rev),
#             peak.RV.Velocity=maxabs(revV),
#             peak.LH.Velocity=maxabs(lev),
#             peak.LV.Velocity=maxabs(levV),
#             peakFR=max(sdf10),
#             post.saccade= first(post.saccade),
#             nspk=sum(rasters),
#             avgFR=nspk/dur*1000,
#             r.amp=sqrt(conj.H.Amp^2+conj.V.Amp^2),
#             asleep=sd.conj.velocity>7.5 || dur>2000,
#             post.saccade=first(post.saccade),
#             burst.index= peak.during.saccade-post.saccade)%>%
#   ungroup() %>%
#   mutate(sacID=row_number())->
#   zsp

t%>% 
  filter(dsnum>0) %>%
  group_by(monkey,neuron,dsnum_extended) %>%
  filter(n()>80,n()<300) %>%
  summarize(dur=n(),
            dsnum=first(dsnum),
            peak.during.saccade=max(sdf10[1:(n()-80)]),
            post.saccade=mean(sdf10[(n()-50):n()]),
            burst.index= peak.during.saccade-post.saccade,
            conj.vert=first(conj.vert),
            upward=conj.vert>0,
            sd.conj.velocity=sd(conj.velocity),
            mean.conj.velocity=mean(conj.velocity),
            sd.verg.velocity=sd(verg.velocity),
            spread=max(conj.velocity)-min(conj.velocity),
            qrange=quantile(conj.velocity,0.975)-quantile(conj.velocity,0.025),
            R.H.Amp=last(rep)-first(rep),
            R.V.Amp=last(repV)-first(repV),
            L.H.Amp=last(lep)-first(lep),
            L.V.Amp=last(lepV)-first(lepV),
            disjunctiveH=sign(R.H.Amp*L.H.Amp)<0,
            disjunctiveV=sign(R.V.Amp*L.V.Amp)<0,
            conj.H.Amp=(R.H.Amp+L.H.Amp)/2,
            conj.V.Amp=(R.V.Amp+L.V.Amp)/2,
            peak.conj.velocity=maxabs(conj.velocity),
            peak.H.Velocity=maxabs((rev+lev)/2),
            peak.V.Velocity=maxabs((revV+levV)/2),
            peak.RH.Velocity=maxabs(rev),
            peak.RV.Velocity=maxabs(revV),
            peak.LH.Velocity=maxabs(lev),
            peak.LV.Velocity=maxabs(levV),
            peakFR=max(sdf10),
            nspk=sum(rasters),
            avgFR=nspk/dur*1000,
            r.amp=sqrt(conj.H.Amp^2+conj.V.Amp^2),
            asleep=sd.conj.velocity>7.5 || dur>2000) ->
  zsp

burst <- filter(zsp,dur<150)

#measure each saccade
burst %>% 
  group_by(neuron,dsnum) %>%
  filter(abs(peak.V.Velocity)<1000,
         abs(peak.H.Velocity)<1000) %>%
  summarize_at(vars(peakFR,
                    peak.V.Velocity,
                    peak.H.Velocity,
                    R.H.Amp,
                    L.H.Amp,
                    R.V.Amp,
                    L.V.Amp,
                    avgFR,
                    disjunctiveH,
                    disjunctiveV,
                    burst.index,
                    post.saccade),
               funs(first)) %>%
  separate(neuron,c('monkey','cellnum'),remove=FALSE)->
  bvelplot

#fit a model for each neuron
burst %>%
  group_by(neuron) %>%
  do(tidy(lm(peakFR~peak.H.Velocity+peak.V.Velocity,data=.))) %>%
  mutate(term=replace(term,term=='(Intercept)','b')) %>%
  select(term,estimate) %>%
  spread(term,estimate) %>%
  separate(neuron,c('monkey','cellnum'),remove=FALSE)->
  bvellm

ver.plot<-bvelplot%>%
  ggplot()+
  geom_point(aes(peak.V.Velocity,peakFR-post.saccade,color=monkey),alpha=0.1)+
  stat_smooth(aes(peak.V.Velocity,peakFR-post.saccade,group=neuron),color='black',method='lm',se = FALSE)+
  facet_wrap(~monkey)+
  ylim(0,NA)+
  xlab('Peak Vertical Eye Velocity (deg/s)')+
  ylab('Firing Rate during Saccade (spk/s)')+
  theme(legend.position = 'none')

hor.plot<-bvelplot%>%
  ggplot()+
  geom_point(aes(peak.H.Velocity,peakFR-post.saccade,color=monkey),alpha=0.1)+
  stat_smooth(aes(peak.H.Velocity,peakFR-post.saccade,group=neuron),color='black',method='lm',se = FALSE)+
  facet_wrap(~monkey)+
  ylim(0,NA)+
  xlab('Peak Horizontal Eye Velocity (deg/s)')+
  ylab('Firing Rate during Saccade (spk/s)')+
  theme(legend.position = 'bottom')

multiplot(ver.plot,hor.plot)

By plotting peak firing rate as a function of peak eye velocity, we can identify cells that have a burst during saccades. These can be distinguished from cells with high tonic firing rates by their steep slope in these plots.

From this plot we can infer that the cells are much more sensitive to vertical eye movements than horizontal eye movements, which is expected for cells in the INC. In addition, none of the cells from DC (the exotrope) show a robust sensitivity to velocity in either the horizontal or vertical directions.

bvellm %>%
  select(neuron,peak.V.Velocity,peak.H.Velocity) %>%
  arrange(desc(abs(peak.V.Velocity))) %>%
  head() %>%
  kable()

neuron	peak.V.Velocity	peak.H.Velocity
Bee-116	-0.2255224	0.0612237
Bee-109	-0.1843909	0.0024562
Kopachuck-908	-0.1807773	-0.0388253
Kopachuck-906	-0.1789044	0.0437980
Bee-110	-0.1707814	0.0205807
Kopachuck-935	0.1445146	0.0100951

Let’s look at some of these cells:

plotsaccades(bvelplot,'Kopachuck-906')

# plotsaccades(bvelplot,'Bee-116')
plotsaccades(bvelplot,'Bee-110')

Looking at the bottom-right plot, there seems to be one cell from Kopachuck (the esotrope) that shows a strong sensitivity to the horizontal component of saccades.

plotsaccades(bvelplot,'Kopachuck-905')

# bvellmplot %>%
#   filter(neuron=='Kopachuck-905') %>%
#   ggplot()+
#   geom_point(aes(peak.H.Velocity,peak.V.Velocity,color=avgFR),size=3)

Let’s compare this will how the cell appeared in our fixation analysis:

multiplot(plotlist=show1cell(zp,zmc,'Kopachuck-905',color1='purple'),cols=2)

This cell has very little sensitivity to eye position during fixation. It appears that it is a single example of a burst neuron. Burst neurons without tonic activity have been reported from this area before, but we haven’t typically encountered many of them. It is also intereting that this cell prefers leftward saccades rather than upward or downard cells, as would be expected for this area. perhaps this is an abnormality related to the strabismus.

#The next step here will be to analyze the preferred direction based on the burst. There should be some effort to distinguish burst-tonic cells from tonic-only or burst-only cells. I'm not sure that we have any burst only cells in this collection, but the algorithm should be able to distinguish them if they should exist. I think I have the code written already to do a lot of this...

#There should be some effort to assess for the possiblity of monocular - bring in some of the analysis I did on the disjunctive saccades. 

#Write up the introduction and methods
measureCellSaccades<-function(x){
  #first I'm just going to make this for the conjugate then I'll come back 
  #and add monocular models if needed

  #Conjugate
  mburst<-lm(peakFR~peak.H.Velocity+peak.V.Velocity,data=x)
  V<-coef(mburst)[['peak.V.Velocity']]
  H<-coef(mburst)[['peak.H.Velocity']]
  
  #calculate direction preference using relative importance of horizontal and vertical
  br=as.numeric(calc.relimp(mburst)$lmg)
  dir.imp<- atan2(br[2]*sign(V),br[1]*sign(H))*180/pi
  
  R2<-summary(mburst)$r.squared
  
  #Right Eye
  mburst<-lm(peakFR~peak.RH.Velocity+peak.RV.Velocity,data=x)
  VR<-coef(mburst)[['peak.RV.Velocity']]
  HR<-coef(mburst)[['peak.RH.Velocity']]
  
  #calculate direction preference using relative importance of horizontal and vertical
  br=as.numeric(calc.relimp(mburst)$lmg)
  dir.imp.R<- atan2(br[2]*sign(VR),br[1]*sign(HR))*180/pi
  
  R2R<-summary(mburst)$r.squared
  
  #Left Eye
  mburst<-lm(peakFR~peak.LH.Velocity+peak.LV.Velocity,data=x)
  VL<-coef(mburst)[['peak.LV.Velocity']]
  HL<-coef(mburst)[['peak.LH.Velocity']]
  
  #calculate direction preference using relative importance of horizontal and vertical
  br=as.numeric(calc.relimp(mburst)$lmg)
  dir.imp.L<- atan2(br[2]*sign(VL),br[1]*sign(HL))*180/pi
  
  R2L<-summary(mburst)$r.squared

  d<-data.frame(neuron=x$neuron[1],
                dir.imp=dir.imp,R2=R2,
                dir.imp.R=dir.imp.R,R2R=R2R,
                dir.imp.L=dir.imp.L,R2L=R2L)
}

measureCellSaccadesBI<-function(x){
  #first I'm just going to make this for the conjugate then I'll come back 
  #and add monocular models if needed

  #Conjugate
  mburst<-lm(burst.index~peak.H.Velocity+peak.V.Velocity,data=x)
  V<-coef(mburst)[['peak.V.Velocity']]
  H<-coef(mburst)[['peak.H.Velocity']]
  
  #calculate direction preference using relative importance of horizontal and vertical
  br=as.numeric(calc.relimp(mburst)$lmg)
  dir.imp<- atan2(br[2]*sign(V),br[1]*sign(H))*180/pi
  
  R2<-summary(mburst)$r.squared
  
  #Right Eye
  mburst<-lm(burst.index~peak.RH.Velocity+peak.RV.Velocity,data=x)
  VR<-coef(mburst)[['peak.RV.Velocity']]
  HR<-coef(mburst)[['peak.RH.Velocity']]
  
  #calculate direction preference using relative importance of horizontal and vertical
  br=as.numeric(calc.relimp(mburst)$lmg)
  dir.imp.R<- atan2(br[2]*sign(VR),br[1]*sign(HR))*180/pi
  
  R2R<-summary(mburst)$r.squared
  
  #Left Eye
  mburst<-lm(burst.index~peak.LH.Velocity+peak.LV.Velocity,data=x)
  VL<-coef(mburst)[['peak.LV.Velocity']]
  HL<-coef(mburst)[['peak.LH.Velocity']]
  
  #calculate direction preference using relative importance of horizontal and vertical
  br=as.numeric(calc.relimp(mburst)$lmg)
  dir.imp.L<- atan2(br[2]*sign(VL),br[1]*sign(HL))*180/pi
  
  R2L<-summary(mburst)$r.squared

  d<-data.frame(neuron=x$neuron[1],
                dir.imp=dir.imp,R2=R2,
                dir.imp.R=dir.imp.R,R2R=R2R,
                dir.imp.L=dir.imp.L,R2L=R2L)
}

burst %>%
  group_by(neuron) %>%
  do(measureCellSaccades(.)) %>%
  separate(neuron,c('monkey','cellnum'),remove=FALSE)->
bsummary

burst %>%
  group_by(neuron) %>%
  do(measureCellSaccadesBI(.)) %>%
  separate(neuron,c('monkey','cellnum'),remove=FALSE)->
bsummaryBI

As above, we calculate the preferred direction of these cells based on the importance of peak horizontal vs vertical conjugate eye velocity. First, we will plot the direction preferences alone, then we scale the lines by the goodness-of-fit, as we did with the fixation analysis.

multiplot(plotlist=list(
ggplot(bsummary)+
  geom_segment(aes(y=0,yend=1,x=dir.imp,xend=dir.imp,color=monkey))+
  scale_x_continuous(limits=c(-180,180),breaks=c(0,90,180,-90))+
  coord_polar(start=pi/2,direction=-1)+
  theme_minimal()+
  theme(legend.position='none')+
  facet_wrap(~monkey,ncol=1)+
  xlab('Direction preference\nusing relative importance')
,
ggplot(bsummary)+
  geom_segment(aes(y=0,yend=R2,x=dir.imp,xend=dir.imp,color=monkey))+
  scale_x_continuous(limits=c(-180,180),breaks=c(0,90,180,-90))+
  coord_polar(start=pi/2,direction=-1)+
  theme_minimal()+
  theme(legend.position='none')+
  facet_wrap(~monkey,ncol=1)+
  xlab('Direction preference\nusing relative importance')
),cols=2)

It appears that the results are very similar to the directions obtained by analysis during fixations. The main exception is that the \(R^2\) values are lower for many cells since only a minority show a distinguishable burst. Notice that the odd strongly horizontal cell from Kopachuck appears in this plot as well, verifying that our methods are reasonable.

Burst Index

The “Burst index” has been used to quantify the strength of the burst component of burst-tonic cells. To calculate the burst index we subtract the tonic firing rate after the saccade from the peak firing rate during the saccade. Cells that do not burst will have a small burst index, while cells that pause during saccades will have a negative burst index. Since cells may burst or pause only for saccades of certain directions or amplitudes, we will have to analyze burst index as a function of saccade metrics.

pvert<-burst %>%
  filter(abs(conj.V.Amp)>5,abs(conj.V.Amp)<40)%>%
  ggplot()+
  stat_smooth(aes(conj.V.Amp,peak.during.saccade-post.saccade,group=neuron),method='lm',se=FALSE)+
  facet_wrap(~monkey)

phor<-burst %>%
  filter(abs(conj.H.Amp)>5,abs(conj.H.Amp)<40)%>%
  ggplot()+
  stat_smooth(aes(conj.H.Amp,peak.during.saccade-post.saccade,group=neuron),method='lm',se=FALSE)+
  facet_wrap(~monkey)

multiplot(pvert,phor,cols=1)

Does using the burst index rather than the peak firing rate change the calculated direction preferences? Burst index is just the peak firing rate minus the post-saccade tonic firing rate.

multiplot(plotlist=list(
ggplot(bsummaryBI)+
  geom_segment(aes(y=0,yend=1,x=dir.imp,xend=dir.imp,color=monkey))+
  scale_x_continuous(limits=c(-180,180),breaks=c(0,90,180,-90))+
  coord_polar(start=pi/2,direction=-1)+
  theme_minimal()+
  theme(legend.position='none')+
  facet_wrap(~monkey,ncol=1)+
  xlab('Direction preference\nusing relative importance')
,
ggplot(bsummaryBI)+
  geom_segment(aes(y=0,yend=R2,x=dir.imp,xend=dir.imp,color=monkey))+
  scale_x_continuous(limits=c(-180,180),breaks=c(0,90,180,-90))+
  coord_polar(start=pi/2,direction=-1)+
  theme_minimal()+
  theme(legend.position='none')+
  facet_wrap(~monkey,ncol=1)+
  xlab('Direction preference\nusing relative importance')
),cols=2)

Monocular analysis

bvelplot%>%
  filter(disjunctiveH) %>%
  ggplot()+
  geom_point(aes(peak.V.Velocity,avgFR,color=monkey),alpha=0.1)+
  stat_smooth(aes(peak.V.Velocity,avgFR,group=neuron),color='black',method='lm',se = FALSE)+
  facet_wrap(~monkey)+
  ylim(0,NA)+
  xlab('Peak Vertical Eye Velocity (deg/s)')+
  ylab('Firing Rate during Saccade (spk/s)')+
  theme(legend.position = 'none')

# bvelplot %>%
#   filter(disjunctiveV,monkey %in% c('Kopachuck','DC')) %>%
#   # filter(disjunctive,monkey=='Kopachuck') %>%
#   group_by(monkey,neuron) %>%
#   do(m=lm(peakFR~R.V.Amp+L.V.Amp,data=.)) %>%
#   mutate(a=as.numeric(calc.relimp(m)$lmg)[1],
#          b=as.numeric(calc.relimp(m)$lmg)[2]) ->
#   bm
# 
# ggplot(bm)+
#   geom_point(aes(a,b,color=monkey))+
#   geom_abline()
# 
# qplot(a-b,data=bm,bins=50)
# 
# filter(bm,abs(a-b)>0.15)

One explanation for the disjunctive eye movements evoked by stimulation is that when we stimulate, we only activate a portion of the neurons in the area. However, if all of the cells are conjugate, then this explanation will work for movements with abnormal directions, but doesn’t explain disconjugate movements. For this explanation to work for disconjugate movements, we would predict that the cells should have different preferences for left or right eyes… Could this also show up as a difference in the weights for the two eyes? Like there are some cells that are strongly right eye and strongly left eye. This would probably show up in the \(R^2\) for the fit.

bsummary %>%
  mutate(r.skew=90-abs(dir.imp.R), 
         l.skew=90-abs(dir.imp.L)) ->
  bsummary

skew.plot.full<- ggplot(bsummary) +
  geom_point(aes(abs(r.skew),abs(l.skew),color=monkey),size=3)+
  geom_abline()+
  ggtitle('Showing cells where left or right eye are skewed less than 20')


skew.plot.zoom<- ggplot(bsummary %>% filter(abs(r.skew)<25 | abs(l.skew)<25)) +
  geom_point(aes(abs(r.skew),abs(l.skew),color=monkey),size=3)+
  geom_abline()

skew.plot.zoom2<- ggplot(bsummary %>% filter(abs(r.skew)<25, abs(l.skew)<25)) +
  geom_point(aes(abs(r.skew),abs(l.skew),color=monkey),size=3)+
  geom_abline()+
  ggtitle('Showing cells where both left and right eye are skewed less than 20')

multiplot(skew.plot.full,skew.plot.zoom,skew.plot.zoom2)

ggplot(bsummary)+
  geom_point(aes(R2R,R2L,color=monkey))+
  geom_abline()

Given how poorly the DC cells fit for the saccadic fits, should we be doing this for the static fixation analysis instead/in addition? Perhaps that will be the next thing to do.

zm %>%
  mutate(r.skew=90-abs(dir.imp.R), 
         l.skew=90-abs(dir.imp.L)) ->
  zm

skew.plot.full<- ggplot(zm) +
  geom_point(aes(abs(r.skew),abs(l.skew),color=monkey),size=3)+
  geom_abline()+
  ggtitle('Showing cells where left or right eye are skewed less than 20')

skew.plot.full

skew.plot<- ggplot(zm) +
  geom_point(aes(r.skew,l.skew,color=monkey),size=3)+
  geom_abline()+
  ggtitle('Showing cells where left or right eye are skewed less than 20')

skew.plot

# ggplotly(skew.plot.full)
# skew.plot.zoom<- ggplot(zm %>% filter(abs(r.skew)<25 | abs(l.skew)<25)) +
#   geom_point(aes(abs(r.skew),abs(l.skew),color=monkey),size=3)+
#   geom_abline()
# 
# skew.plot.zoom2<- ggplot(zm %>% filter(abs(r.skew)<25, abs(l.skew)<25)) +
#   geom_point(aes(abs(r.skew),abs(l.skew),color=monkey),size=3)+
#   geom_abline()+
#   ggtitle('Showing cells where both left and right eye are skewed less than 20')
# 
# multiplot(skew.plot.full,skew.plot.zoom,skew.plot.zoom2)

ggplot(zm)+
  geom_point(aes(R2R,R2L,color=monkey))+
  geom_abline()

Analysis of INC in strabismus

Adam Pallus