rmeddis@38: function testRP2(MAPparamsName,paramChanges)
rmeddis@38: % testIHC evaluates IHC I/O function
rmeddis@38: % multiple BF can be used but only one is easier to interpret.
rmeddis@38: % e.g. testRP(1000,'Normal',{});
rmeddis@38: 
rmeddis@38: global experiment method inputStimulusParams
rmeddis@38: global stimulusParameters IHC_VResp_VivoParams IHC_cilia_RPParams
rmeddis@38: savePath=path;
rmeddis@38: addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])
rmeddis@38: dbstop if error
rmeddis@38: 
rmeddis@38: figure(4), clf,
rmeddis@38: set (gcf, 'name', ['IHC'])
rmeddis@38: set(gcf,'position',[613   354   360   322])
rmeddis@38: drawColors='rgbkmcy';
rmeddis@38: drawnow
rmeddis@38: 
rmeddis@38: if nargin<3
rmeddis@38:     paramChanges=[];
rmeddis@38: end
rmeddis@38: if nargin<2
rmeddis@38:     MAPparamsName='Normal';
rmeddis@38: end
rmeddis@38: if nargin<3
rmeddis@38:     BF=800;
rmeddis@38: end
rmeddis@38: 
rmeddis@38: levels=-20:10:100;
rmeddis@38: nLevels=length(levels);
rmeddis@38: toneDuration=.05;
rmeddis@38: silenceDuration=.01;
rmeddis@38: sampleRate=50000;
rmeddis@38: dt=1/sampleRate;
rmeddis@38: 
rmeddis@38: allIHC_RP_peak=[];
rmeddis@38: allIHC_RP_dc=[];
rmeddis@38: 
rmeddis@38: %% Ruggero
rmeddis@38: %%Ruggero data
rmeddis@38: RuggeroData=[
rmeddis@38:     0	2.00E-10;
rmeddis@38:     10	5.00E-10;
rmeddis@38:     20	1.50E-09;
rmeddis@38:     30	2.50E-09;
rmeddis@38:     40	5.30E-09;
rmeddis@38:     50	1.00E-08;
rmeddis@38:     60	1.70E-08;
rmeddis@38:     70	2.50E-08;
rmeddis@38:     80	4.00E-08;
rmeddis@38:     90	6.00E-08;
rmeddis@38:     100	1.50E-07;
rmeddis@38:     110	3.00E-07;
rmeddis@38:     ];
rmeddis@38: 
rmeddis@38: 
rmeddis@38: BF=10000;
rmeddis@38: targetFrequency=BF;
rmeddis@38: 
rmeddis@38: IHC_RP_peak=zeros(nLevels,1);
rmeddis@38: IHC_RP_min=zeros(nLevels,1);
rmeddis@38: IHC_RP_dc=zeros(nLevels,1);
rmeddis@38: 
rmeddis@38: time=dt:dt:toneDuration;
rmeddis@38: 
rmeddis@38: rampDuration=0.004;
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: toneStartptr=length(silence)+1;
rmeddis@38: toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
rmeddis@38: toneEndptr=toneStartptr+round(toneDuration/dt) -1;
rmeddis@38: 
rmeddis@38: levelNo=0;
rmeddis@38: for leveldB=levels
rmeddis@38:     levelNo=levelNo+1;
rmeddis@38:     % replicate at all levels
rmeddis@38:     amp=28e-6*10^(leveldB/20);
rmeddis@38:     
rmeddis@38:     %% create signal (leveldB/ frequency)
rmeddis@38:     inputSignal=amp*sin(2*pi*targetFrequency'*time);
rmeddis@38:     inputSignal= ramp.*inputSignal;
rmeddis@38:     inputSignal=[silence inputSignal silence];
rmeddis@38:     inputStimulusParams.sampleRate=1/dt;
rmeddis@38:     %         global IHC_ciliaParams
rmeddis@38:     
rmeddis@38:     %% disable efferent for fast processing
rmeddis@38:     
rmeddis@38:     %% run the model
rmeddis@38:     global  DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
rmeddis@38:         IHCoutput
rmeddis@38:     AN_spikesOrProbability='probability';
rmeddis@38:     
rmeddis@38:     MAP1_14(inputSignal, sampleRate, BF, ...
rmeddis@38:         MAPparamsName, AN_spikesOrProbability, paramChanges);
rmeddis@38:     
rmeddis@38:     % DRNL
rmeddis@38:     DRNLoutput=DRNLoutput;
rmeddis@38:     DRNL_peak(levelNo,1)=max(DRNLoutput(toneMidptr:toneEndptr));
rmeddis@38:     DRNL_min(levelNo,1)=min(DRNLoutput(toneMidptr:toneEndptr));
rmeddis@38:     DRNL_dc(levelNo,1)=mean(DRNLoutput(toneMidptr:toneEndptr));
rmeddis@38:     
rmeddis@38: end
rmeddis@38: %% plot DRNL
rmeddis@38: subplot(2,2,1)
rmeddis@38: referenceDisp=10e-9;
rmeddis@38: semilogy(levels,DRNL_peak, 'linewidth',2), hold on
rmeddis@38: semilogy(RuggeroData(:,1),RuggeroData(:,2),'o')
rmeddis@38: title(['BM: Ruggero ' num2str(BF) ' Hz'])
rmeddis@38: ylabel ('displacement(m)'), xlabel('dB SPL')
rmeddis@38: xlim([min(levels) max(levels)]), ylim([1e-10 1e-7])
rmeddis@38: grid on
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% Dallos
rmeddis@38: BF=800;
rmeddis@38: targetFrequency=BF;
rmeddis@38: 
rmeddis@38: IHC_RP_peak=zeros(nLevels,1);
rmeddis@38: IHC_RP_min=zeros(nLevels,1);
rmeddis@38: IHC_RP_dc=zeros(nLevels,1);
rmeddis@38: 
rmeddis@38: time=dt:dt:toneDuration;
rmeddis@38: 
rmeddis@38: rampDuration=0.004;
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: toneStartptr=length(silence)+1;
rmeddis@38: toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
rmeddis@38: toneEndptr=toneStartptr+round(toneDuration/dt) -1;
rmeddis@38: 
rmeddis@38: levelNo=0;
rmeddis@38: for leveldB=levels
rmeddis@38:     levelNo=levelNo+1;
rmeddis@38:     % replicate at all levels
rmeddis@38:     amp=28e-6*10^(leveldB/20);
rmeddis@38:     
rmeddis@38:     %% create signal (leveldB/ frequency)
rmeddis@38:     inputSignal=amp*sin(2*pi*targetFrequency'*time);
rmeddis@38:     inputSignal= ramp.*inputSignal;
rmeddis@38:     inputSignal=[silence inputSignal silence];
rmeddis@38:     inputStimulusParams.sampleRate=1/dt;
rmeddis@38:     %         global IHC_ciliaParams
rmeddis@38:     
rmeddis@38:     %% disable efferent for fast processing
rmeddis@38:     method.DRNLSave=1;
rmeddis@38:     method.IHC_cilia_RPSave=1;
rmeddis@38:     method.IHCpreSynapseSave=1;
rmeddis@38:     method.IHC_cilia_RPSave=1;
rmeddis@38:     method.segmentDuration=-1;
rmeddis@38:     moduleSequence=1:4;
rmeddis@38:     
rmeddis@38:     %% run the model
rmeddis@38:     global  DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
rmeddis@38:         IHCoutput
rmeddis@38:     AN_spikesOrProbability='probability';
rmeddis@38:     
rmeddis@38:     MAP1_14(inputSignal, sampleRate, BF, ...
rmeddis@38:         MAPparamsName, AN_spikesOrProbability, paramChanges);
rmeddis@38:     
rmeddis@38:     % DRNL
rmeddis@38:     DRNLoutput=DRNLoutput;
rmeddis@38:     DRNL_peak(levelNo,1)=max(DRNLoutput(toneMidptr:toneEndptr));
rmeddis@38:     DRNL_min(levelNo,1)=min(DRNLoutput(toneMidptr:toneEndptr));
rmeddis@38:     DRNL_dc(levelNo,1)=mean(DRNLoutput(toneMidptr:toneEndptr));
rmeddis@38:     
rmeddis@38:     % cilia
rmeddis@38:     IHC_ciliaData=IHC_cilia_output;
rmeddis@38:     IHC_ciliaData=IHC_ciliaData;
rmeddis@38:     IHC_cilia_peak(levelNo,1)=...
rmeddis@38:         max(IHC_ciliaData(toneMidptr:toneEndptr));
rmeddis@38:     IHC_cilia_min(levelNo,1)=...
rmeddis@38:         min(IHC_ciliaData(toneMidptr:toneEndptr));
rmeddis@38:     IHC_cilia_dc(levelNo,1)=...
rmeddis@38:         mean(IHC_ciliaData(toneMidptr:toneEndptr));
rmeddis@38:     
rmeddis@38:     % RP
rmeddis@38:     IHC_RPData=IHCoutput;
rmeddis@38:     IHC_RPData=IHC_RPData;
rmeddis@38:     IHC_RP_peak(levelNo,1)=...
rmeddis@38:         max(IHC_RPData(toneMidptr:toneEndptr));
rmeddis@38:     IHC_RP_min(levelNo,1)=...
rmeddis@38:         min(IHC_RPData(toneMidptr:toneEndptr));
rmeddis@38:     IHC_RP_dc(levelNo,1)=...
rmeddis@38:         mean(IHC_RPData(toneMidptr:toneEndptr));
rmeddis@38:     
rmeddis@38: end
rmeddis@38: 
rmeddis@38: %% Dallos and Harris data
rmeddis@38: %% plot receptor potentials
rmeddis@38: figure(4)
rmeddis@38: subplot(2,2,3)
rmeddis@38: % RP I/O function min and max
rmeddis@38: restingRP=IHC_RP_peak(1);
rmeddis@38: toPlot= [fliplr(IHC_RP_min(:,1)') IHC_RP_peak(:,1)'];
rmeddis@38: microPa=   28e-6*10.^(levels/20);
rmeddis@38: microPa=[-fliplr(microPa) microPa];
rmeddis@38: plot(microPa,toPlot, 'linewidth',2)
rmeddis@38: 
rmeddis@38: dallosx=[-0.9	-0.1	-0.001	0.001	0.01	0.9];
rmeddis@38: dallosy=[-8	-7.8	-6.5	11	16.5	22]/1000 + restingRP;
rmeddis@38: subplot(2,2,3)
rmeddis@38: hold on, plot(dallosx,dallosy, 'o')
rmeddis@38: plot([-1 1], [restingRP restingRP], 'r')
rmeddis@38: title(' Dallos(86) data at 800 Hz')
rmeddis@38: ylabel ('receptor potential(V)'), xlabel('Pa')
rmeddis@38: ylim([-0.08 -0.02]), xlim([-1 1])
rmeddis@38: grid on
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% Patuzzi
rmeddis@38: BF=7000;
rmeddis@38: targetFrequency=BF;
rmeddis@38: 
rmeddis@38: IHC_RP_peak=zeros(nLevels,1);
rmeddis@38: IHC_RP_min=zeros(nLevels,1);
rmeddis@38: IHC_RP_dc=zeros(nLevels,1);
rmeddis@38: 
rmeddis@38: time=dt:dt:toneDuration;
rmeddis@38: 
rmeddis@38: rampDuration=0.004;
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: toneStartptr=length(silence)+1;
rmeddis@38: toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
rmeddis@38: toneEndptr=toneStartptr+round(toneDuration/dt) -1;
rmeddis@38: 
rmeddis@38: levelNo=0;
rmeddis@38: for leveldB=levels
rmeddis@38:     levelNo=levelNo+1;
rmeddis@38:     % replicate at all levels
rmeddis@38:     amp=28e-6*10^(leveldB/20);
rmeddis@38:     
rmeddis@38:     %% create signal (leveldB/ frequency)
rmeddis@38:     inputSignal=amp*sin(2*pi*targetFrequency'*time);
rmeddis@38:     inputSignal= ramp.*inputSignal;
rmeddis@38:     inputSignal=[silence inputSignal silence];
rmeddis@38:     inputStimulusParams.sampleRate=1/dt;
rmeddis@38:     %         global IHC_ciliaParams
rmeddis@38:     
rmeddis@38:     %% run the model
rmeddis@38:     global  DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
rmeddis@38:         IHCoutput
rmeddis@38:     AN_spikesOrProbability='probability';
rmeddis@38:     
rmeddis@38:     MAP1_14(inputSignal, sampleRate, BF, ...
rmeddis@38:         MAPparamsName, AN_spikesOrProbability, paramChanges);
rmeddis@38:     
rmeddis@38:     % RP
rmeddis@38:     IHC_RPData=IHCoutput;
rmeddis@38:     IHC_RP_peak(levelNo,1)=...
rmeddis@38:         max(IHC_RPData(toneMidptr:toneEndptr));
rmeddis@38:     IHC_RP_min(levelNo,1)=...
rmeddis@38:         min(IHC_RPData(toneMidptr:toneEndptr));
rmeddis@38:     IHC_RP_dc(levelNo,1)=...
rmeddis@38:         mean(IHC_RPData(toneMidptr:toneEndptr));
rmeddis@38: end
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% RP I/O function min and max
rmeddis@38: figure(4)
rmeddis@38: subplot(2,2,4)
rmeddis@38: restingRP=IHC_RP_peak(1);
rmeddis@38: peakRP=max(IHC_RP_peak);
rmeddis@38: plot(levels, IHC_RP_peak, 'linewidth',2)
rmeddis@38: hold on
rmeddis@38: plot(levels, IHC_RP_dc, ':', 'linewidth',2)
rmeddis@38: plot([min(levels) max(levels)], [restingRP restingRP], 'r')
rmeddis@38: xlim([min(levels) max(levels)])
rmeddis@38: % animal data
rmeddis@38: sndLevel=[5	15	25	35	45	55	65	75];
rmeddis@38: RPanimal=restingRP+[0.5	2	4.6	5.8	6.4	7.2	8	10.2]/1000;
rmeddis@38: % could be misleading when restingRP changes
rmeddis@38: RPanimal=-0.060+[0.5	2	4.6	5.8	6.4	7.2	8	10.2]/1000;
rmeddis@38: hold on, plot(sndLevel,RPanimal,'o')
rmeddis@38: 
rmeddis@38: grid on
rmeddis@38: title(' 7 kHz Patuzzi')
rmeddis@38: ylabel ('RP(V) peak and DC'), xlabel('dB SPL')
rmeddis@38: ylim([-0.07 -0.04])
rmeddis@38: path(savePath);
rmeddis@38: disp(paramChanges)