rmeddis@38: function vectorStrength=testAN(probeFrequency,BFlist, levels, ...
rmeddis@38:     paramsName,paramChanges)
rmeddis@38: % testAN generates rate/level functions for AN and brainstem units.
rmeddis@38: %  also other information like PSTHs, MOC efferent activity levels.
rmeddis@38: % Vector strength calculations require the computation of period
rmeddis@38: % histograms. for this reason the sample rate must always be an integer
rmeddis@38: % multiple of the tone frequency. This applies to both the sampleRate and
rmeddis@38: % the spikes sampling rate.
rmeddis@38: % A full 'spikes' model is used.
rmeddis@38: % paramChanges is a cell array of strings containing MATLAB statements that
rmeddis@38: % change model parameters. See MAPparamsNormal for examples.
rmeddis@38: % e.g.
rmeddis@38: % testAN(1000,1000, -10:10:80,'Normal',{});
rmeddis@35: 
rmeddis@38: 
rmeddis@38: global IHC_VResp_VivoParams  IHC_cilia_RPParams IHCpreSynapseParams
rmeddis@38: global AN_IHCsynapseParams
rmeddis@38: global ANoutput dtSpikes CNoutput ICoutput ICmembraneOutput ANtauCas
rmeddis@38: global ARattenuation MOCattenuation
rmeddis@38: tic
rmeddis@36: dbstop if error
rmeddis@36: restorePath=path;
rmeddis@38: addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ...
rmeddis@38:     ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...
rmeddis@38:     ['..' filesep 'testPrograms'])
rmeddis@35: 
rmeddis@38: if nargin<5, paramChanges={'IHCpreSynapseParams.tauCa= [30e-6 90e-6];'}; end
rmeddis@38: if nargin<4, paramsName='PL'; end
rmeddis@38: if nargin<3, levels=-10:10:80; end
rmeddis@38: if nargin==0, probeFrequency=1000; BFlist=1000; end
rmeddis@38: nLevels=length(levels);
rmeddis@36: 
rmeddis@38: toneDuration=.2;   rampDuration=0.005;   silenceDuration=.02;
rmeddis@38: localPSTHbinwidth=0.001;
rmeddis@36: 
rmeddis@38: %% guarantee that the sample rate is an interger multiple 
rmeddis@38: %   of the probe frequency
rmeddis@38: % we want 5 bins per period for spikes
rmeddis@38: spikesSampleRate=5*probeFrequency;
rmeddis@38: % model sample rate must be an integer multiple of this and in the region
rmeddis@38: % of 50000
rmeddis@38: sampleRate=spikesSampleRate*round(50000/spikesSampleRate);
rmeddis@38: dt=1/sampleRate;
rmeddis@38: % avoid very slow spikes sampling rate
rmeddis@38: spikesSampleRate=spikesSampleRate*ceil(10000/spikesSampleRate);
rmeddis@36: 
rmeddis@38: %% pre-allocate storage
rmeddis@38: AN_HSRonset=zeros(nLevels,1);
rmeddis@38: AN_HSRsaturated=zeros(nLevels,1);
rmeddis@38: AN_LSRonset=zeros(nLevels,1);
rmeddis@38: AN_LSRsaturated=zeros(nLevels,1);
rmeddis@38: CNLSRrate=zeros(nLevels,1);
rmeddis@38: CNHSRsaturated=zeros(nLevels,1);
rmeddis@38: ICHSRsaturated=zeros(nLevels,1);
rmeddis@38: ICLSRsaturated=zeros(nLevels,1);
rmeddis@38: vectorStrength=zeros(nLevels,1);
rmeddis@36: 
rmeddis@38: AR=zeros(nLevels,1);
rmeddis@38: MOC=zeros(nLevels,1);
rmeddis@36: 
rmeddis@38: figure(15), clf
rmeddis@38: set(gcf,'position',[980   356   401   321])
rmeddis@38: figure(5), clf
rmeddis@38: set(gcf,'position', [980 34 400 295])
rmeddis@38: drawnow
rmeddis@36: 
rmeddis@36: 
rmeddis@38: %% main computational loop (vary level)
rmeddis@38: levelNo=0;
rmeddis@38: for leveldB=levels
rmeddis@38:     levelNo=levelNo+1;
rmeddis@38:     amp=28e-6*10^(leveldB/20);
rmeddis@38:     fprintf('%4.0f\t', leveldB)
rmeddis@36: 
rmeddis@38:     %% generate tone and silences
rmeddis@38:     time=dt:dt:toneDuration;
rmeddis@38:     rampTime=dt:dt:rampDuration;
rmeddis@38:     ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
rmeddis@38:         ones(1,length(time)-length(rampTime))];
rmeddis@38:     ramp=ramp.*fliplr(ramp);
rmeddis@38:     
rmeddis@38:     silence=zeros(1,round(silenceDuration/dt));
rmeddis@38:     
rmeddis@38:     inputSignal=amp*sin(2*pi*probeFrequency'*time);
rmeddis@38:     inputSignal= ramp.*inputSignal;
rmeddis@38:     inputSignal=[silence inputSignal];
rmeddis@38:     
rmeddis@38:     %% run the model
rmeddis@38:     AN_spikesOrProbability='spikes';  
rmeddis@38:     nExistingParamChanges=length(paramChanges);
rmeddis@38:     paramChanges{nExistingParamChanges+1}=...
rmeddis@38:         ['AN_IHCsynapseParams.spikesTargetSampleRate=' ...
rmeddis@38:         num2str(spikesSampleRate) ';'];
rmeddis@36: 
rmeddis@38:     MAP1_14(inputSignal, 1/dt, BFlist, ...
rmeddis@38:         paramsName, AN_spikesOrProbability, paramChanges);
rmeddis@38:     
rmeddis@38:     nTaus=length(ANtauCas);
rmeddis@38:     
rmeddis@38:     %% Auditory nerve evaluate and display (Fig. 5)
rmeddis@38:     %LSR (same as HSR if no LSR fibers present)
rmeddis@38:     [nANFibers nTimePoints]=size(ANoutput);
rmeddis@38:     numLSRfibers=nANFibers/nTaus;
rmeddis@38:     numHSRfibers=numLSRfibers;
rmeddis@38:     
rmeddis@38:     LSRspikes=ANoutput(1:numLSRfibers,:);
rmeddis@38:     PSTH=UTIL_PSTHmaker(LSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38:     PSTHLSR=mean(PSTH,1)/localPSTHbinwidth;  % across fibers rates
rmeddis@38:     PSTHtime=localPSTHbinwidth:localPSTHbinwidth:...
rmeddis@38:         localPSTHbinwidth*length(PSTH);
rmeddis@38:     AN_LSRonset(levelNo)= max(PSTHLSR); % peak in 5 ms window
rmeddis@38:     AN_LSRsaturated(levelNo)= mean(PSTHLSR(round(length(PSTH)/2):end));
rmeddis@38:     
rmeddis@38:     % AN HSR
rmeddis@38:     HSRspikes= ANoutput(end- numHSRfibers+1:end, :);
rmeddis@38:     PSTH=UTIL_PSTHmaker(HSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38:     PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)
rmeddis@38:     AN_HSRonset(levelNo)= max(PSTH);
rmeddis@38:     AN_HSRsaturated(levelNo)= mean(PSTH(round(length(PSTH)/2): end));
rmeddis@38:     
rmeddis@38:     figure(5), subplot(2,2,2)
rmeddis@38:     hold off, bar(PSTHtime,PSTH, 'b')
rmeddis@38:     hold on,  bar(PSTHtime,PSTHLSR,'r')
rmeddis@38:     ylim([0 1000])
rmeddis@38:     xlim([0 length(PSTH)*localPSTHbinwidth])
rmeddis@38:     grid on
rmeddis@38:     set(gcf,'name',['PSTH: ' num2str(BFlist), ' Hz: ' num2str(leveldB) ' dB']);
rmeddis@38:     
rmeddis@38:     % AN - CV
rmeddis@38:     %  CV is computed 5 times. Use the middle one (3) as most typical
rmeddis@38:     cvANHSR=  UTIL_CV(HSRspikes, dtSpikes);
rmeddis@38:     
rmeddis@38:     % AN - vector strength
rmeddis@38:     PSTH=sum(CNoutput (11:20,:));
rmeddis@38:     [PH, binTimes]=UTIL_periodHistogram...
rmeddis@38:         (PSTH, dtSpikes, probeFrequency);
rmeddis@38:     VS=UTIL_vectorStrength(PH);
rmeddis@38:     vectorStrength(levelNo)=VS;
rmeddis@38:     disp(['sat rate= ' num2str(AN_HSRsaturated(levelNo)) ...
rmeddis@38:         ';   phase-locking VS = ' num2str(VS)])
rmeddis@38:     title(['AN HSR: CV=' num2str(cvANHSR(3),'%5.2f') ...
rmeddis@38:         'VS=' num2str(VS,'%5.2f')])
rmeddis@38:     
rmeddis@38:     %% CN - first-order neurons
rmeddis@38:     
rmeddis@38:     % CN driven by LSR fibers
rmeddis@38:     [nCNneurons c]=size(CNoutput);
rmeddis@38:     nLSRneurons=round(nCNneurons/nTaus);
rmeddis@38:     CNLSRspikes=CNoutput(1:nLSRneurons,:);
rmeddis@38:     PSTH=UTIL_PSTHmaker(CNLSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38:     PSTH=sum(PSTH)/nLSRneurons;
rmeddis@38: %     CNLSRrate(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth;
rmeddis@38:     CNLSRrate(levelNo)=mean(PSTH)/localPSTHbinwidth;
rmeddis@38:     
rmeddis@38:     %CN HSR
rmeddis@38:     MacGregorMultiHSRspikes=...
rmeddis@38:         CNoutput(end-nLSRneurons+1:end,:);
rmeddis@38:     PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38:     PSTH=sum(PSTH)/nLSRneurons;
rmeddis@38:     PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)
rmeddis@38:     
rmeddis@38: %     CNHSRsaturated(levelNo)=mean(PSTH(length(PSTH)/2:end));
rmeddis@38:     CNHSRsaturated(levelNo)=mean(PSTH);
rmeddis@38:     
rmeddis@38:     figure(5), subplot(2,2,3)
rmeddis@38:     bar(PSTHtime,PSTH)
rmeddis@38:     ylim([0 1000])
rmeddis@38:     xlim([0 length(PSTH)*localPSTHbinwidth])
rmeddis@38:     cvMMHSR= UTIL_CV(MacGregorMultiHSRspikes, dtSpikes);
rmeddis@38:     title(['CN    CV= ' num2str(cvMMHSR(3),'%5.2f')])
rmeddis@38:     
rmeddis@38:     %% IC LSR
rmeddis@38:     [nICneurons c]=size(ICoutput);
rmeddis@38:     nLSRneurons=round(nICneurons/nTaus);
rmeddis@38:     ICLSRspikes=ICoutput(1:nLSRneurons,:);
rmeddis@38:     PSTH=UTIL_PSTHmaker(ICLSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38: %     ICLSRsaturated(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth;
rmeddis@38:     ICLSRsaturated(levelNo)=mean(PSTH)/localPSTHbinwidth;
rmeddis@38:     
rmeddis@38:     %IC HSR
rmeddis@38:     MacGregorMultiHSRspikes=...
rmeddis@38:         ICoutput(end-nLSRneurons+1:end,:);
rmeddis@38:     PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38:     ICHSRsaturated(levelNo)= (sum(PSTH)/nLSRneurons)/toneDuration;
rmeddis@36: 
rmeddis@38:     AR(levelNo)=min(ARattenuation);
rmeddis@38:     MOC(levelNo)=min(MOCattenuation(length(MOCattenuation)/2:end));
rmeddis@38:     
rmeddis@38:     time=dt:dt:dt*size(ICmembraneOutput,2);
rmeddis@38:     figure(5), subplot(2,2,4)
rmeddis@38:     % plot HSR (last of two)
rmeddis@38:     plot(time,ICmembraneOutput(end-nLSRneurons+1, 1:end),'k')
rmeddis@38:     ylim([-0.07 0])
rmeddis@38:     xlim([0 max(time)])
rmeddis@38:     title(['IC  ' num2str(leveldB,'%4.0f') 'dB'])
rmeddis@38:     drawnow
rmeddis@38:     
rmeddis@38:     figure(5), subplot(2,2,1)
rmeddis@38:     plot(20*log10(MOC), 'k'), hold on
rmeddis@38:     plot(20*log10(AR), 'r'), hold off
rmeddis@38:     title(' MOC'), ylabel('dB attenuation')
rmeddis@38:     ylim([-30 0])
rmeddis@38:     
rmeddis@38: %     x=input('prompt')
rmeddis@38: end % level
rmeddis@36: 
rmeddis@38: %% plot with levels on x-axis
rmeddis@38: figure(5), subplot(2,2,1)
rmeddis@38: plot(levels,20*log10(MOC), 'k'),hold on
rmeddis@38: plot(levels,20*log10(AR), 'r'), hold off
rmeddis@38: title(' MOC'), ylabel('dB attenuation')
rmeddis@38: ylim([-30 0])
rmeddis@38: xlim([0 max(levels)])
rmeddis@36: 
rmeddis@38: fprintf('\n')
rmeddis@38: toneDuration=2;
rmeddis@38: rampDuration=0.004;
rmeddis@38: silenceDuration=.02;
rmeddis@38: nRepeats=200;   % no. of AN fibers
rmeddis@38: fprintf('toneDuration  %6.3f\n', toneDuration)
rmeddis@38: fprintf(' %6.0f  AN fibers (repeats)\n', nRepeats)
rmeddis@38: fprintf('levels')
rmeddis@38: fprintf('%6.2f\t', levels)
rmeddis@38: fprintf('\n')
rmeddis@36: 
rmeddis@38: %% ---------------------- display summary results (Fig 15)
rmeddis@38: figure(15), clf
rmeddis@38: nRows=2; nCols=2;
rmeddis@36: 
rmeddis@38: % AN rate - level ONSET functions
rmeddis@38: subplot(nRows,nCols,1)
rmeddis@38: plot(levels,AN_LSRonset,'ro'), hold on
rmeddis@38: plot(levels,AN_HSRonset,'ko'), hold off
rmeddis@38: ylim([0 1000])
rmeddis@38: if length(levels)>1
rmeddis@38:     xlim([min(levels) max(levels)])
rmeddis@36: end
rmeddis@36: 
rmeddis@38: ttl=['tauCa= ' num2str(IHCpreSynapseParams.tauCa)];
rmeddis@38: title( ttl)
rmeddis@38: xlabel('level dB SPL'), ylabel('peak rate (sp/s)'), grid on
rmeddis@38: text(0, 800, 'AN onset', 'fontsize', 14)
rmeddis@38: 
rmeddis@38: % AN rate - level ADAPTED function
rmeddis@38: subplot(nRows,nCols,2)
rmeddis@38: plot(levels,AN_LSRsaturated, 'ro'), hold on
rmeddis@38: plot(levels,AN_HSRsaturated, 'ko'), hold off
rmeddis@38: ylim([0 400])
rmeddis@38: set(gca,'ytick',0:50:300)
rmeddis@38: if length(levels)>1
rmeddis@38: xlim([min(levels) max(levels)])
rmeddis@36: end
rmeddis@38: set(gca,'xtick',[levels(1):20:levels(end)])
rmeddis@38: %     grid on
rmeddis@38: ttl=[   'spont=' num2str(mean(AN_HSRsaturated(1,:)),'%4.0f')...
rmeddis@38:     '  sat=' num2str(mean(AN_HSRsaturated(end,1)),'%4.0f')];
rmeddis@38: title( ttl)
rmeddis@38: xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
rmeddis@38: text(0, 340, 'AN adapted', 'fontsize', 14), grid on
rmeddis@36: 
rmeddis@38: % CN rate - level function
rmeddis@38: subplot(nRows,nCols,3)
rmeddis@38: plot(levels,CNLSRrate, 'ro'), hold on
rmeddis@38: plot(levels,CNHSRsaturated, 'ko'), hold off
rmeddis@38: ylim([0 400])
rmeddis@38: set(gca,'ytick',0:50:300)
rmeddis@38: if length(levels)>1
rmeddis@38: xlim([min(levels) max(levels)])
rmeddis@38: end
rmeddis@38: set(gca,'xtick',[levels(1):20:levels(end)])
rmeddis@38: %     grid on
rmeddis@38: ttl=[   'spont=' num2str(mean(CNHSRsaturated(1,:)),'%4.0f') '  sat=' ...
rmeddis@38:     num2str(mean(CNHSRsaturated(end,1)),'%4.0f')];
rmeddis@38: title( ttl)
rmeddis@38: xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
rmeddis@38: text(0, 350, 'CN', 'fontsize', 14), grid on
rmeddis@36: 
rmeddis@38: % IC rate - level function
rmeddis@38: subplot(nRows,nCols,4)
rmeddis@38: plot(levels,ICLSRsaturated, 'ro'), hold on
rmeddis@38: plot(levels,ICHSRsaturated, 'ko'), hold off
rmeddis@38: ylim([0 400])
rmeddis@38: set(gca,'ytick',0:50:300)
rmeddis@38: if length(levels)>1
rmeddis@38: xlim([min(levels) max(levels)])
rmeddis@38: end
rmeddis@38: set(gca,'xtick',[levels(1):20:levels(end)]), grid on
rmeddis@38: ttl=['spont=' num2str(mean(ICHSRsaturated(1,:)),'%4.0f') ...
rmeddis@38:     '  sat=' num2str(mean(ICHSRsaturated(end,1)),'%4.0f')];
rmeddis@38: title( ttl)
rmeddis@38: xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
rmeddis@38: text(0, 350, 'IC', 'fontsize', 14)
rmeddis@38: set(gcf,'name',' AN CN IC rate/level')
rmeddis@36: 
rmeddis@38: fprintf('\n')
rmeddis@38: disp('levels vectorStrength')
rmeddis@38: fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength'])
rmeddis@38: fprintf('\n')
rmeddis@38: fprintf('Phase locking, max vector strength=\t %6.4f\n\n',...
rmeddis@38:     max(vectorStrength))
rmeddis@36: 
rmeddis@38: allData=[ levels'  AN_HSRonset AN_HSRsaturated...
rmeddis@38:     AN_LSRonset AN_LSRsaturated ...
rmeddis@38:     CNHSRsaturated CNLSRrate...
rmeddis@38:     ICHSRsaturated ICLSRsaturated];
rmeddis@38: fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR  \tICLSR \n');
rmeddis@38: UTIL_printTabTable(round(allData))
rmeddis@38: fprintf('VS (phase locking)= \t%6.4f\n\n',...
rmeddis@38:     max(vectorStrength))
rmeddis@36: 
rmeddis@38: UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams')
rmeddis@38: UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')
rmeddis@38: UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')
rmeddis@38: 
rmeddis@35: fprintf('\n')
rmeddis@38: disp('levels vectorStrength')
rmeddis@38: fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength'])
rmeddis@38: fprintf('\n')
rmeddis@38: fprintf('Phase locking, max vector strength= \t%6.4f\n\n',...
rmeddis@38:     max(vectorStrength))
rmeddis@35: 
rmeddis@38: allData=[ levels'  AN_HSRonset AN_HSRsaturated...
rmeddis@38:     AN_LSRonset AN_LSRsaturated ...
rmeddis@38:     CNHSRsaturated CNLSRrate...
rmeddis@38:     ICHSRsaturated ICLSRsaturated];
rmeddis@38: fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR  \tICLSR \n');
rmeddis@38: UTIL_printTabTable(round(allData))
rmeddis@38: fprintf('VS (phase locking)= \t%6.4f\n\n',...
rmeddis@38:     max(vectorStrength))
rmeddis@36: 
rmeddis@38: path(restorePath)
rmeddis@36: toc