Spatio-temporal LN model.ipf

#pragma rtGlobals=1		// Use modern global access method.

//Function CalcRFfn
//Calculates receptive field by correlating a spike train with a stimulus
//This is the spike-triggered average, the stimulus preceding a spike is summed,
//then divided by the number of spikes. If you want to calculate the best fit
//linear model (linear kernel), and interpret the amplitude of the receptive field,
//you must multiply afterwards by the firing rate. If you are calculating an LN model,
//multiplying by the firing rate is not necessary.
//Note, this function currently does not subtract the mean stimulus (an option in Correlate)
//This only becomes necessary when mapping the receptive field using membrane potential
//artificially encoded as a spike train to work with Correlate. The very large number of spikes
//created with this method (hundreds of Hz) makes it necessary to subtract this offset.
//Input variables:
//st:Wave containing white noise stimulus
//sp:Wave containing spike train
//tstart:Start time of receptive field relative to a spike
//tend:End time of receptive field relative to a spike
//tstep:Time step for receptive field
//Output:
//wave containing receptive field is '"name of spike train" + "_rf"
function CalcRFfn(st,sp,tstart,tend,tstep)
	wave st, sp
	variable tstart, tend, tstep
	variable tau, i, pnt
	//Wave to hold receptive field
	make /o/n=(dimsize(st,0),dimsize(st,1),(tend-tstart)/tstep) $nameofwave(sp)+"_rf"
	wave rf=$nameofwave(sp)+"_rf"
	setscale /p z,tstart,tstep,rf
	make /o/n=(dimsize(st,0),dimsize(st,1)) onefr //One frame of receptive field
	variable lastframe=dimoffset(st,2)+dimdelta(st,2)*dimsize(st,2)
	for (tau=tstart;tau<tend;tau+=tstep)  //Loop over time delay and calculate
		onefr=0						   //spike-triggered average stimulus
		for (i=0;(i<numpnts (sp)) %& (sp[i]<(lastframe-tend));i+=1)
			onefr+=st[p][q](sp[i]+tau)
		endfor
		rf[][][pnt]=onefr[p][q]
		pnt+=1
	endfor
	rf/=numpnts(sp)					 //Divide by number of spikes
	wavestats /q rf;rf-=v_avg			 //Subtract off mean value of entire RF
	killwaves /z onefr
end

//Macro STLNmodel
//This function creates a spatio-temporal LN model
//INPUT VARIABLES:
//sstim:White noise stimulus
//resp:Response of cell
//rfin:Linear spatio-temporal receptive field of cell
macro SpatiotemporalLNmodel(sstim,sresp,srfin)
	string sstim,sresp,srfin
	Prompt sstim,"White noise stimulus"
	Prompt sresp,"Spiking response"
	Prompt srfin,"Receptive field"
	variable pnum
	variable threshval=3	// Value for setting receptive field threshold
	variable tstart=-0.025 //Start and end of time interval used to judge
	variable tend = -0.125  //which pixels are significant
	threshrf ($srfin,threshval,tstart,tend) //Threshold receptive field to take only significant pixels
	convolverf ($sstim,$srfin,rfx,rfy) //Convolve rf with stim
	wavestats /q convsum 
	//Convert spike train from resp into a function of time by binning in 1 ms intervals
	spikehist(sresp,rightx(convsum),0.001,0,1,"no","no")
	//Place linear prediction (in wave 'convsum') into wave 'lp' and response into wave 'v'
	variable npoints=min(numpnts(convsum),numpnts($sresp+"_h")) //Make sure lp and v have
	duplicate /o/r=[0,npoints-1] $sresp+"_h",v						//same number of points
	duplicate /o/r=[0,npoints-1] convsum,lp		
	//Calculate nonlinearity with 10000 (or numpnts(lp)) points
	dosnl (min(10000,numpnts(lp))) 
	duplicate /o vpy,$srfin+".snl1y"
	duplicate /o vpx,$srfin+".snl1x"
	duplicate /o vpmean $srfin+".snl1"
	dowindow /f binsnlgraph //Bring SNL binning window to front
end

//Function threshrf
//Threshold receptive field
//INPUT VARIABLES
//rfin: wave containing receptive field
//thr:threshold value in standard deviations of the overall receptive field value
//tstart:starting time of interval used to judge which pixels are significant
//tend: ending time of interval used  to judge which pixels are significant
function threshrf (rfin,thr,tstart,tend)
	wave rfin
	variable thr
	variable tstart
	variable tend
	//rfmags:root-mean-squared value  of each pixel
	//rfx:wave containing x coord of significant pixels
	//rfy:wave containing y coord of significant pixels
	make /o/n=(dimsize (rfin,0)*dimsize(rfin,1)) rfmags,rfx,rfy
	variable px,py,pix
	make /o/n=(dimsize(rfin,2)) onepixtimesq // For one pixel p, holds p(t)^2, used to
	                                                                         // compute the rms value of that pixel over time
	setscale /p x,dimoffset(rfin,2),dimdelta(rfin,2),onepixtimesq
	for (px=0;px<dimsize(rfin,0);px+=1)        //Loop over x coord
		for (py=0;py<dimsize(rfin,1);py+=1) //Loop over y coord
			onepixtimesq=rfin[px][py][p]^2
			wavestats /q/r=(tstart,tend) onepixtimesq 
			rfmags[pix]=sqrt(v_avg) //rms value of one pix
			rfx[pix]=px
			rfy[pix]=py
			pix+=1
		endfor
	endfor
	sort /r rfmags,rfmags,rfx,rfy //sort the pixels with the greatest rms values
	wavestats /q  rfmags
	rfx=selectnumber( rfmags[p]>v_sdev*thr,nan,rfx) //get pixels greater than the thr in std devs
	rfy=selectnumber( rfmags[p]>v_sdev*thr,nan,rfy) //get pixels greater than the thr in std devs
	rfmags=selectnumber( rfmags[p]>v_sdev*thr,nan,rfmags)
	wavestats /q rfmags
	deletepoints v_npnts,numpnts(rfmags),rfmags,rfx,rfy 
	make /o/n=(dimsize(rfin,0),dimsize(rfin,1)) rf_thr //thresholded spatial receptive field
	rf_thr=0
	for (pix=0;pix<v_npnts;pix+=1)
		rf_thr[rfx[pix]][rfy[pix]]=rfmags[pix]
	endfor
end

//function convolverf
//Convolves the significant pixels of a spatio-temporal receptive field with a stimulus
//Each pixel of the receptive field is convolved with the corresponding pixel of the stimulus
//resulting in a single function of time for each pixel. All single pixel convolutions
//are then summed yielding a single function of time for the linear prediction of a response.
//INPUT VARIABLES
//wst:wave containing stimulus
//wrf:wave containing receptive field
//wrfx:wave containing list of x coordinates of significant pixels
//wrfy:wave containing list of y coordinates of significant pixels
function convolverf (wst,wrf,wrfx,wrfy)
	wave wst
	wave wrf
	wave wrfx
	wave wrfy
	variable pix
	//st1p:stimulus one pixel in milliseconds, number of points is
	//(number of frames in stimulus) * (frame time in seconds) * (1000 ms / s)
	//cor1p: convolution for one pixel, points are in ms
	//corsum: sum of all pixel 
	make /o/n=(dimsize(wst,2)) st1p,conv1p,convsum
	setscale /p x,0,dimdelta(wst,2),st1p,conv1p,convsum
	//rf1p:time course of one pixel of receptive field in ms, number of points is
	//(number of timepoints in receptive field) * (delta-t in seconds) * (1000 ms / s)
	make /o/n=(dimsize(wrf,2)) rf1p
	setscale /p x,0,dimdelta(wrf,2),rf1p
	convsum=0
	for (pix=0;pix<numpnts (wrfx);pix+=1) //Loop over number of significant pixels in RF
		rf1p=wrf[wrfx[pix]][wrfy[pix]](-x)
		st1p=wst[wrfx[pix]][wrfy[pix]](x)
		duplicate /o st1p,conv1p
		convolve rf1p,conv1p
		convsum+=conv1p
		deletepoints 0,numpnts(rf1p),conv1p
	endfor
end

macro spikehist(Nwavein,periodlength,binlength,starttime,numperiods,onflag,displayflag)
	string Nwavein,onflag="no",displayflag="no"
	variable periodlength=60,binlength=0.6,starttime=60,numperiods=34
	prompt Nwavein,"",popup,Wavelist("*_s",";","")+""
	prompt onflag,"On and off responses?",popup,"yes;no"
	prompt displayflag,"Display result?",popup,"yes;no"
	PauseUpdate; Silent 1		// building window...
	duplicate /o $Nwavein spikesbin
	copylimit(Nwavein,"spikesbin",starttime,periodlength*numperiods+starttime)
	//spikesbin/=1000
	spikesbin = mod(spikesbin,periodlength)
	string Nwaveinpref=stringfromlist(0,Nwavein,"_")
	make /o/n=(periodlength/binlength) histwave
	histwave=0
	Histogram/B={0,binlength,(periodlength/binlength)} spikesbin, histwave
	setscale /p x,binlength/2,binlength,histwave
	histwave/=numperiods*binlength
	wavestats /q histwave
	if (stringmatch(onflag,"yes"))
		string Nwaveouton,Nwaveoutoff
		sprintf Nwaveouton,  Nwaveinpref+"_onh"
		sprintf Nwaveoutoff,  Nwaveinpref+"_offh"	
		make /o/n=(numpnts(histwave)/2) $Nwaveouton,$Nwaveoutoff
		$Nwaveouton=0
		$Nwaveoutoff=0
		setscale /p x,binlength/2,binlength*2,$Nwaveouton
		setscale /p x,binlength*3/2,binlength*2,$Nwaveoutoff
		variable i=0
		do
			$Nwaveouton[i]=histwave[i*2]
			$Nwaveoutoff[i]=histwave[i*2+1]
			i+=1
		while (i<numpnts($Nwaveouton))
		//dowindow /k gspikehist
		display $Nwaveouton,$Nwaveoutoff as Nwaveinpref+ " spike histogram"
		//dowindow /c gspikehist	
		ModifyGraph mode=3
		ModifyGraph marker=19
		SetAxis left -0.1,v_max*1.1 
		ModifyGraph rgb($Nwaveouton)=(0,0,65535)
	else
		string Nwaveout
		sprintf Nwaveout,  Nwavein+"_h"
		duplicate /o histwave,$Nwaveout
	endif	
	if (stringmatch(displayflag,"yes"))
		
		display $Nwaveout as Nwaveinpref+ " spike histogram"
		ModifyGraph mode=3
		ModifyGraph marker=19
		Textbox/N=text0/F=0/A=MC "\\F'Helvetica'\\Z14"+Nwavein
		Textbox/C/N=text0/X=-40.00/Y=44.81
	endif
endmacro

function DoSNL(numbins)
	variable numbins
	silent 1;pauseupdate
	variable binpoints,stepfactor=2
	numbins/=stepfactor
	duplicate /o lp,lpwave // Linear prediction
	duplicate /o v,Vwave   // Response
	sort lpwave,lpwave,vwave
	binpoints=numpnts(vwave)/numbins
	make /o/n=(numpnts(vwave)) vpn //SNLmean,num in each bin
	make /o/n=(numpnts(vwave)) vpmean,Vpsd,vpx,vpy,bst,bend // SNL sd
	vpn=0
	vpsd=0
	wavestats /q lpwave
	setscale /i x,v_min,v_max,vpmean,vpn,vpsd,vpx,vpy
	variable bin,binstart,binend
	binstart=0
	bin=1
	duplicate /o vwave,cfit
	cfit=0
	do	
		binend=binstart+binpoints-1
		//Make sure the last bin isn't too small
		if ((numpnts(lpwave)-binend)<binpoints/stepfactor)
			binend=numpnts(lpwave)
		endif
		bst[bin]=binstart
		bend[bin]=binend
		wavestats /q/r=[binstart,binend] vwave
		vpy[bin]=v_avg
		vpsd[bin]=v_sdev
		vpn[bin]=binpoints
		wavestats /q/r=[binstart,binend] lpwave
		vpx[bin]=v_avg
		//First bin
		if (bin==1)
			bst[bin]=binstart
			bend[bin]=binend
			vpy[0]=cfit[binstart]
			vpsd[0]=v_sdev
			vpn[0]=binpoints
			vpx[0]=lpwave[binstart]
		endif
		bin=bin+1
		binstart=binstart+binpoints/stepfactor
	while (binstart<=(numpnts(lpwave)-binpoints))
	duplicate /o vpsd,vpsem
	vpsem=vpsd/sqrt(vpn)
	deletepoints bin,numpnts(vwave),vpy,vpsd,vpx,vpsem
	//Interpolate SNL into evenly spaced points
	make /o/n=1000 snlmean
	wavestats /q vpx
	setscale /i x,v_min,v_max,snlmean
	execute "interpolate /y=snlmean vpy /x=vpx"
	//killwaves /z vwave,lpwave,vpx,vpy
end