Mercurial > hg > map
changeset 26:b03ef38fe497
refractory (prob) effect added
| author | Ray Meddis <rmeddis@essex.ac.uk> |
|---|---|
| date | Tue, 21 Jun 2011 14:58:12 +0100 |
| parents | d2c4c07df02c |
| children | d4a7675b0413 |
| files | Help and reference data/MAP1_14 quick reference.doc Help and reference data/MAPlog.doc MAP/MAP1_14.m MAP/MAP1_14parallel.m multithreshold 1.46/old files/MAPmodel.m multithreshold 1.46/savedData/mostRecentResults.mat multithreshold 1.46/testAN.m multithreshold 1.46/testBM.m multithreshold 1.46/testSynapse.m parameterStore/MAPparamsEndo.m parameterStore/MAPparamsNormal.m parameterStore/MAPparamsNormalTest.m parameterStore/MAPparamsOHC.m parameterStore/MAPparamsStd.m testPrograms/demoTwisterProbability.m testPrograms/demoTwisterSpikes.m testPrograms/test_MAP1_14.m utilities/UTIL_PSTHmaker.m utilities/UTIL_showMAP.m |
| diffstat | 19 files changed, 511 insertions(+), 145 deletions(-) [+] |
line wrap: on
line diff
--- a/MAP/MAP1_14.m Fri Jun 17 14:30:28 2011 +0100 +++ b/MAP/MAP1_14.m Tue Jun 21 14:58:12 2011 +0100 @@ -52,29 +52,26 @@ % AR is computed using across channel activity % MOC is computed on a within-channel basis. +if nargin<1 + error(' MAP1_14 is not a script but a function that must be called') +end + +if nargin<6 + paramChanges=[]; +end +% Read parameters from MAPparams<***> file in 'parameterStore' folder +cmd=['method=MAPparams' MAPparamsName ... + '(BFlist, sampleRate, 0, paramChanges);']; +eval(cmd); +% Beware, 'BFlist=-1' is a legitimate argument for MAPparams<> +% if the calling program allows MAPparams to specify the list +BFlist=DRNLParams.nonlinCFs; % save as global for later plotting if required savedBFlist=BFlist; saveAN_spikesOrProbability=AN_spikesOrProbability; saveMAPparamsName=MAPparamsName; -% Read parameters from MAPparams<***> file in 'parameterStore' folder -cmd=['method=MAPparams' MAPparamsName ... - '(BFlist, sampleRate, 0);']; -eval(cmd); - -% Beware, 'BFlist=-1' is a legitimate argument for MAPparams<> -% if the calling program allows MAPparams to specify the list -BFlist=DRNLParams.nonlinCFs; - -% now accept last mintue parameter changes required by the calling program -if nargin>5 && ~isempty(paramChanges) - nChanges=length(paramChanges); - for idx=1:nChanges - eval(paramChanges{idx}) - end -end - dt=1/sampleRate; duration=length(inputSignal)/sampleRate; % segmentDuration is specified in parameter file (must be >efferent delay) @@ -190,7 +187,7 @@ MOCrateToAttenuationFactor=DRNLParams.rateToAttenuationFactor; rateToAttenuationFactorProb=DRNLParams.rateToAttenuationFactorProb; -MOCrateThreshold=DRNLParams.MOCrateThreshold; +MOCrateThresholdProb=DRNLParams.MOCrateThresholdProb; % smoothing filter for MOC a1=dt/DRNLParams.MOCtau-1; a0=1; @@ -355,6 +352,7 @@ % Spikes are necessary if CN and IC are to be computed nFibersPerChannel= AN_IHCsynapseParams.numFibers; nANfibers= nChannels*nFibersPerChannel; +AN_refractory_period= AN_IHCsynapseParams.refractory_period; y=AN_IHCsynapseParams.y; l=AN_IHCsynapseParams.l; @@ -377,12 +375,12 @@ ANprobability=zeros(nChannels,segmentLength); ANprobRateOutput=zeros(nChannels,signalLength); +lengthAbsRefractoryP= round(AN_refractory_period/dt); % special variables for monitoring synaptic cleft (specialists only) savePavailableSeg=zeros(nChannels,segmentLength); savePavailable=zeros(nChannels,signalLength); % spikes % ! ! ! ! ! ! ! ! -AN_refractory_period= AN_IHCsynapseParams.refractory_period; lengthAbsRefractory= round(AN_refractory_period/ANdt); AN_ydt= repmat(AN_IHCsynapseParams.y*ANdt, nANfibers,1); @@ -777,7 +775,7 @@ [smoothedRates, MOCprobBoundary{idx}] = ... filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ... MOCprobBoundary{idx}); - smoothedRates=smoothedRates-MOCrateThreshold; + smoothedRates=smoothedRates-MOCrateThresholdProb; smoothedRates(smoothedRates<0)=0; MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ... (1- smoothedRates* rateToAttenuationFactorProb); @@ -1131,4 +1129,24 @@ end % segment +%% apply refractory correction to spike probabilities + +% switch AN_spikesOrProbability +% case 'probability' +% ANprobOutput=ANprobRateOutput*dt; +% [r c]=size(ANprobOutput); +% % find probability of no spikes in refractory period +% pNoSpikesInRefrac=ones(size(ANprobOutput)); +% pSpike=zeros(size(ANprobOutput)); +% for epochNo=lengthAbsRefractoryP+2:c +% pNoSpikesInRefrac(:,epochNo)=... +% pNoSpikesInRefrac(:,epochNo-2)... +% .*(1-pSpike(:,epochNo-1))... +% ./(1-pSpike(:,epochNo-lengthAbsRefractoryP-1)); +% pSpike(:,epochNo)= ANprobOutput(:,epochNo)... +% .*pNoSpikesInRefrac(:,epochNo); +% end +% ANprobRateOutput=pSpike/dt; +% end + path(restorePath)
--- a/MAP/MAP1_14parallel.m Fri Jun 17 14:30:28 2011 +0100 +++ b/MAP/MAP1_14parallel.m Tue Jun 21 14:58:12 2011 +0100 @@ -192,7 +192,7 @@ MOCrateToAttenuationFactor=DRNLParams.rateToAttenuationFactor; rateToAttenuationFactorProb=DRNLParams.rateToAttenuationFactorProb; -MOCrateThreshold=DRNLParams.MOCrateThreshold; +MOCrateThresholdProb=DRNLParams.MOCrateThresholdProb; % smoothing filter for MOC % Nyquist=(1/ANdt)/2; @@ -819,7 +819,7 @@ [smoothedRates, MOCprobBoundary{idx}] = ... filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ... MOCprobBoundary{idx}); - smoothedRates=smoothedRates-MOCrateThreshold; + smoothedRates=smoothedRates-MOCrateThresholdProb; smoothedRates(smoothedRates<0)=0; MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ... (1- smoothedRates* rateToAttenuationFactorProb);
--- a/multithreshold 1.46/old files/MAPmodel.m Fri Jun 17 14:30:28 2011 +0100 +++ b/multithreshold 1.46/old files/MAPmodel.m Tue Jun 21 14:58:12 2011 +0100 @@ -30,7 +30,7 @@ options.printFiringRates=1; options.showACF=0; options.showEfferent=1; - UTIL_showMAP(options) + UTIL_showMAP(options, paramChanges) end % No response, probably caused by hitting 'stop' button
--- a/multithreshold 1.46/testAN.m Fri Jun 17 14:30:28 2011 +0100 +++ b/multithreshold 1.46/testAN.m Tue Jun 21 14:58:12 2011 +0100 @@ -108,7 +108,7 @@ numHSRfibers=numLSRfibers; LSRspikes=ANoutput(1:numLSRfibers,:); - PSTH=UTIL_makePSTH(LSRspikes, ANdt, localPSTHbinwidth); + PSTH=UTIL_PSTHmaker(LSRspikes, ANdt, localPSTHbinwidth); PSTHLSR=mean(PSTH,1)/localPSTHbinwidth; % across fibers rates PSTHtime=localPSTHbinwidth:localPSTHbinwidth:... localPSTHbinwidth*length(PSTH); @@ -117,7 +117,7 @@ % HSR HSRspikes= ANoutput(end- numHSRfibers+1:end, :); - PSTH=UTIL_makePSTH(HSRspikes, ANdt, localPSTHbinwidth); + PSTH=UTIL_PSTHmaker(HSRspikes, ANdt, localPSTHbinwidth); PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only) AN_HSRonset(levelNo)= max(PSTH); AN_HSRsaturated(levelNo)= mean(PSTH(round(length(PSTH)/2): end)); @@ -150,14 +150,14 @@ [nCNneurons c]=size(CNoutput); nLSRneurons=round(nCNneurons/nTaus); CNLSRspikes=CNoutput(1:nLSRneurons,:); - PSTH=UTIL_makePSTH(CNLSRspikes, ANdt, localPSTHbinwidth); + PSTH=UTIL_PSTHmaker(CNLSRspikes, ANdt, localPSTHbinwidth); PSTH=sum(PSTH)/nLSRneurons; CNLSRrate(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth; %CN HSR MacGregorMultiHSRspikes=... CNoutput(end-nLSRneurons:end,:); - PSTH=UTIL_makePSTH(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth); + PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth); PSTH=sum(PSTH)/nLSRneurons; PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only) @@ -174,13 +174,13 @@ [nICneurons c]=size(ICoutput); nLSRneurons=round(nICneurons/nTaus); ICLSRspikes=ICoutput(1:nLSRneurons,:); - PSTH=UTIL_makePSTH(ICLSRspikes, ANdt, localPSTHbinwidth); + PSTH=UTIL_PSTHmaker(ICLSRspikes, ANdt, localPSTHbinwidth); ICLSRsaturated(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth; %IC HSR MacGregorMultiHSRspikes=... ICoutput(end-nLSRneurons:end,:); - PSTH=UTIL_makePSTH(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth); + PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth); PSTH=sum(PSTH)/nLSRneurons; PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)
--- a/multithreshold 1.46/testBM.m Fri Jun 17 14:30:28 2011 +0100 +++ b/multithreshold 1.46/testBM.m Tue Jun 21 14:58:12 2011 +0100 @@ -19,7 +19,7 @@ % levels= 50; nLevels=length(levels); relativeFrequencies=[0.25 .5 .75 1 1.25 1.5 2]; -% relativeFrequencies=1; +relativeFrequencies=1; % refBMdisplacement is the displacement of the BM at threshold % 1 nm disp at threshold (9 kHz, Ruggero)
--- a/multithreshold 1.46/testSynapse.m Fri Jun 17 14:30:28 2011 +0100 +++ b/multithreshold 1.46/testSynapse.m Tue Jun 21 14:58:12 2011 +0100 @@ -26,6 +26,7 @@ useEfferent=0; % no efferent silenceDuration=0.005; +silenceDuration=0.015; maskerDuration=0.1; recoveryDuration=0.15; rampDuration=0.004;
--- a/parameterStore/MAPparamsEndo.m Fri Jun 17 14:30:28 2011 +0100 +++ b/parameterStore/MAPparamsEndo.m Tue Jun 21 14:58:12 2011 +0100 @@ -1,5 +1,5 @@ function method=MAPparamsEndo ... - (BFlist, sampleRate, showParams) + (BFlist, sampleRate, showParams, paramChanges) % MAPparams<> establishes a complete set of MAP parameters % Parameter file names must be of the form <MAPparams> <name> % @@ -58,11 +58,11 @@ % i.e. a minimum ratio of 0.056. % 'spikes' model: AR based on brainstem spiking activity (LSR) OMEParams.rateToAttenuationFactor=0.006; % * N(all ICspikes) -% OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) +% OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) % 'probability model': Ar based on AN firing probabilities (LSR) OMEParams.rateToAttenuationFactorProb=0.01;% * N(all ANrates) -% OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates) +% OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates) % asymptote should be around 100-200 ms OMEParams.ARtau=.05; % AR smoothing function @@ -104,9 +104,9 @@ DRNLParams.rateToAttenuationFactor = .01; % strength of MOC % DRNLParams.rateToAttenuationFactor = 0; % strength of MOC % 'probability' model: MOC based on AN spiking activity (HSR) -DRNLParams.rateToAttenuationFactorProb = .005; % strength of MOC +DRNLParams.rateToAttenuationFactorProb = .0055; % strength of MOC % DRNLParams.rateToAttenuationFactorProb = .0; % strength of MOC -DRNLParams.MOCrateThreshold =70; % spikes/s probability only +DRNLParams.MOCrateThresholdProb =70; % spikes/s probability only DRNLParams.MOCtau =.1; % smoothing for MOC @@ -263,6 +263,16 @@ end +%% now accept last minute parameter changes required by the calling program +% paramChanges +if nargin>3 && ~isempty(paramChanges) + nChanges=length(paramChanges); + for idx=1:nChanges + eval(paramChanges{idx}) + end +end + + %% write all parameters to the command window % showParams is currently set at the top of htis function if showParams @@ -276,6 +286,14 @@ eval(['UTIL_showStructureSummary(' nm{i} ', ''' nm{i} ''', 10)']) end end + + % highlight parameter changes made locally + if nargin>3 && ~isempty(paramChanges) + fprintf('\n Local parameter changes:\n') + for i=1:length(paramChanges) + disp(paramChanges{i}) + end + end end
--- a/parameterStore/MAPparamsNormal.m Fri Jun 17 14:30:28 2011 +0100 +++ b/parameterStore/MAPparamsNormal.m Tue Jun 21 14:58:12 2011 +0100 @@ -1,5 +1,5 @@ function method=MAPparamsNormal ... - (BFlist, sampleRate, showParams) + (BFlist, sampleRate, showParams, paramChanges) % MAPparams<> establishes a complete set of MAP parameters % Parameter file names must be of the form <MAPparams> <name> % @@ -11,11 +11,11 @@ % output argument % method passes a miscelleny of values -global inputStimulusParams OMEParams DRNLParams +global inputStimulusParams OMEParams DRNLParams IHC_cilia_RPParams global IHC_VResp_VivoParams IHCpreSynapseParams AN_IHCsynapseParams global MacGregorParams MacGregorMultiParams filteredSACFParams global experiment % used by calls from multiThreshold only -global IHC_cilia_RPParams + currentFile=mfilename; % i.e. the name of this mfile method.parameterSource=currentFile(10:end); % for the record @@ -58,11 +58,11 @@ % i.e. a minimum ratio of 0.056. % 'spikes' model: AR based on brainstem spiking activity (LSR) OMEParams.rateToAttenuationFactor=0.006; % * N(all ICspikes) -% OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) +% OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) % 'probability model': Ar based on AN firing probabilities (LSR) OMEParams.rateToAttenuationFactorProb=0.01;% * N(all ANrates) -% OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates) +% OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates) % asymptote should be around 100-200 ms OMEParams.ARtau=.05; % AR smoothing function @@ -104,9 +104,9 @@ DRNLParams.rateToAttenuationFactor = .01; % strength of MOC % DRNLParams.rateToAttenuationFactor = 0; % strength of MOC % 'probability' model: MOC based on AN spiking activity (HSR) -DRNLParams.rateToAttenuationFactorProb = .005; % strength of MOC +DRNLParams.rateToAttenuationFactorProb = .0055; % strength of MOC % DRNLParams.rateToAttenuationFactorProb = .0; % strength of MOC -DRNLParams.MOCrateThreshold =70; % spikes/s probability only +DRNLParams.MOCrateThresholdProb =70; % spikes/s probability only DRNLParams.MOCtau =.1; % smoothing for MOC @@ -263,6 +263,16 @@ end +%% now accept last minute parameter changes required by the calling program +% paramChanges +if nargin>3 && ~isempty(paramChanges) + nChanges=length(paramChanges); + for idx=1:nChanges + eval(paramChanges{idx}) + end +end + + %% write all parameters to the command window % showParams is currently set at the top of htis function if showParams @@ -276,63 +286,15 @@ eval(['UTIL_showStructureSummary(' nm{i} ', ''' nm{i} ''', 10)']) end end -end - - -% ********************************************************************** comparison data -% store individual data here for display on the multiThreshold GUI (if used) -% the final value in each vector is an identifier (BF or duration)) -if isstruct(experiment) - switch experiment.paradigm - case {'IFMC','IFMC_8ms'} - % based on MPa - comparisonData=[ - 66 51 49 48 46 45 54 250; - 60 54 46 42 39 49 65 500; - 64 51 38 32 33 59 75 1000; - 59 51 36 30 41 81 93 2000; - 71 63 53 44 36 76 95 4000; - 70 64 43 35 35 66 88 6000; - 110 110 110 110 110 110 110 8000; - ]; - if length(BFlist)==1 && ~isempty(comparisonData) - availableFrequencies=comparisonData(:,end)'; - findRow= find(BFlist==availableFrequencies); - if ~isempty (findRow) - experiment.comparisonData=comparisonData(findRow,:); - end - end - - case {'TMC','TMC_8ms'} - % based on MPa - comparisonData=[ - 48 58 63 68 75 80 85 92 99 250; - 33 39 40 49 52 61 64 77 79 500; - 39 42 50 81 83 92 96 97 110 1000; - 24 26 32 37 46 51 59 71 78 2000; - 65 68 77 85 91 93 110 110 110 4000; - 20 19 26 44 80 95 96 110 110 6000; - ]; - if length(BFlist)==1 && ~isempty(comparisonData) - availableFrequencies=comparisonData(:,end)'; - findRow= find(BFlist==availableFrequencies); - if ~isempty (findRow) - experiment.comparisonData=comparisonData(findRow,:); - end - end - - case { 'absThreshold', 'absThreshold_8'} - % MPa thresholds - experiment.comparisonData=[ - 32 26 16 18 22 22 0.008; - 16 13 6 9 15 11 0.500 - ]; - - - otherwise - experiment.comparisonData=[]; + % highlight parameter changes made locally + if nargin>3 && ~isempty(paramChanges) + fprintf('\n Local parameter changes:\n') + for i=1:length(paramChanges) + disp(paramChanges{i}) + end end end - +% for backward compatibility +experiment.comparisonData=[];
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/parameterStore/MAPparamsNormalTest.m Tue Jun 21 14:58:12 2011 +0100 @@ -0,0 +1,340 @@ +function method=MAPparamsNormalTest ... + (BFlist, sampleRate, showParams) +% MAPparams<> establishes a complete set of MAP parameters +% Parameter file names must be of the form <MAPparams> <name> +% +% input arguments +% BFlist (optional) specifies the desired list of channel BFs +% otherwise defaults set below +% sampleRate (optional), default is 50000. +% showParams (optional) =1 prints out the complete set of parameters +% output argument +% method passes a miscelleny of values + +global inputStimulusParams OMEParams DRNLParams +global IHC_VResp_VivoParams IHCpreSynapseParams AN_IHCsynapseParams +global MacGregorParams MacGregorMultiParams filteredSACFParams +global experiment % used by calls from multiThreshold only +global IHC_cilia_RPParams + +currentFile=mfilename; % i.e. the name of this mfile +method.parameterSource=currentFile(10:end); % for the record + +efferentDelay=0.010; +method.segmentDuration=efferentDelay; + +if nargin<3, showParams=0; end +if nargin<2, sampleRate=50000; end +if nargin<1 || BFlist(1)<0 % if BFlist= -1, set BFlist to default + lowestBF=250; highestBF= 8000; numChannels=21; + % 21 chs (250-8k)includes BFs at 250 500 1000 2000 4000 8000 + BFlist=round(logspace(log10(lowestBF),log10(highestBF),numChannels)); +end +% BFlist=1000; + +% preserve for backward campatibility +method.nonlinCF=BFlist; +method.dt=1/sampleRate; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% set model parameters +%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% #1 inputStimulus +inputStimulusParams=[]; +inputStimulusParams.sampleRate= sampleRate; + +%% #2 outerMiddleEar +OMEParams=[]; % clear the structure first +% outer ear resonances band pass filter [gain lp order hp] +OMEParams.externalResonanceFilters= [ 10 1 1000 4000]; + +% highpass stapes filter +% Huber gives 2e-9 m at 80 dB and 1 kHz (2e-13 at 0 dB SPL) +OMEParams.OMEstapesLPcutoff= 1000; +OMEParams.stapesScalar= 45e-9; + +% Acoustic reflex: maximum attenuation should be around 25 dB Price (1966) +% i.e. a minimum ratio of 0.056. +% 'spikes' model: AR based on brainstem spiking activity (LSR) +OMEParams.rateToAttenuationFactor=0.006; % * N(all ICspikes) +% OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) + +% 'probability model': Ar based on AN firing probabilities (LSR) +OMEParams.rateToAttenuationFactorProb=0.01;% * N(all ANrates) +% OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates) + +% asymptote should be around 100-200 ms +OMEParams.ARtau=.05; % AR smoothing function +% delay must be longer than the segment length +OMEParams.ARdelay=efferentDelay; %Moss gives 8.5 ms latency +OMEParams.ARrateThreshold=0; + +%% #3 DRNL +DRNLParams=[]; % clear the structure first +DRNLParams.BFlist=BFlist; + +% DRNL nonlinear path +DRNLParams.a=5e4; % DRNL.a=0 means no OHCs (no nonlinear path) +DRNLParams.a=0e4; % DRNL.a=0 means no OHCs (no nonlinear path) + +DRNLParams.b=8e-6; % *compression threshold raised compression +% DRNLParams.b=12e-6; % *compression threshold raised compression +% DRNLParams.b=1; % b=1 means no compression + +DRNLParams.c=0.2; % compression exponent +% nonlinear filters +DRNLParams.nonlinCFs=BFlist; +DRNLParams.nonlinOrder= 3; % order of nonlinear gammatone filters +p=0.2895; q=170; % human (% p=0.14; q=366; % cat) +DRNLParams.nlBWs= p * BFlist + q; +DRNLParams.p=p; DRNLParams.q=q; % save p and q for printing only + +% DRNL linear path: +DRNLParams.g=1000; % linear path gain factor +% linCF is not necessarily the same as nonlinCF +minLinCF=153.13; coeffLinCF=0.7341; % linCF>nonlinBF for BF < 1 kHz +DRNLParams.linCFs=minLinCF+coeffLinCF*BFlist; +DRNLParams.linOrder= 3; % order of linear gammatone filters +minLinBW=100; coeffLinBW=0.6531; +DRNLParams.linBWs=minLinBW + coeffLinBW*BFlist; % bandwidths of linear filters + +% DRNL MOC efferents +DRNLParams.MOCdelay = efferentDelay; % must be < segment length! + +% 'spikes' model: MOC based on brainstem spiking activity (HSR) +DRNLParams.rateToAttenuationFactor = .01; % strength of MOC +% DRNLParams.rateToAttenuationFactor = 0; % strength of MOC +% 'probability' model: MOC based on AN spiking activity (HSR) +DRNLParams.rateToAttenuationFactorProb = .005; % strength of MOC +% DRNLParams.rateToAttenuationFactorProb = .0; % strength of MOC +DRNLParams.MOCrateThresholdProb =70; % spikes/s probability only + +DRNLParams.MOCtau =.1; % smoothing for MOC + + +%% #4 IHC_cilia_RPParams + +IHC_cilia_RPParams.tc= 0.0003; % 0.0003 filter time simulates viscocity +% IHC_cilia_RPParams.tc= 0.0005; % 0.0003 filter time simulates viscocity +IHC_cilia_RPParams.C= 0.01; % 0.1 scalar (C_cilia ) +IHC_cilia_RPParams.u0= 5e-10; +IHC_cilia_RPParams.s0= 10e-10; +IHC_cilia_RPParams.u1= 5e-10; +IHC_cilia_RPParams.s1= 5e-10; + +IHC_cilia_RPParams.Gmax= 4e-9; % 2.5e-9 maximum conductance (Siemens) +IHC_cilia_RPParams.Ga= 1.5e-9; % 4.3e-9 fixed apical membrane conductance + +% #5 IHC_RP +IHC_cilia_RPParams.Cab= 4e-012; % IHC capacitance (F) +IHC_cilia_RPParams.Cab= 2e-012; % IHC capacitance (F) +IHC_cilia_RPParams.Et= 0.100; % endocochlear potential (V) + +IHC_cilia_RPParams.Gk= 2e-008; % 1e-8 potassium conductance (S) +IHC_cilia_RPParams.Ek= -0.08; % -0.084 K equilibrium potential +IHC_cilia_RPParams.Rpc= 0.04; % combined resistances + + +%% #5 IHCpreSynapse +IHCpreSynapseParams=[]; +IHCpreSynapseParams.GmaxCa= 14e-9;% maximum calcium conductance +IHCpreSynapseParams.GmaxCa= 12e-9;% maximum calcium conductance +IHCpreSynapseParams.ECa= 0.066; % calcium equilibrium potential +IHCpreSynapseParams.beta= 400; % determine Ca channel opening +IHCpreSynapseParams.gamma= 100; % determine Ca channel opening +IHCpreSynapseParams.tauM= 0.00005; % membrane time constant ?0.1ms +IHCpreSynapseParams.power= 3; +% reminder: changing z has a strong effect on HF thresholds (like Et) +IHCpreSynapseParams.z= 1e42; % scalar Ca -> vesicle release rate + +LSRtauCa=35e-6; HSRtauCa=85e-6; % seconds +% LSRtauCa=35e-6; HSRtauCa=70e-6; % seconds +IHCpreSynapseParams.tauCa= [LSRtauCa HSRtauCa]; %LSR and HSR fiber + +%% #6 AN_IHCsynapse +% c=kym/(y(l+r)+kl) (spontaneous rate) +% c=(approx) ym/l (saturated rate) +AN_IHCsynapseParams=[]; % clear the structure first +AN_IHCsynapseParams.M= 12; % maximum vesicles at synapse +AN_IHCsynapseParams.y= 4; % depleted vesicle replacement rate +% AN_IHCsynapseParams.y= 6; % depleted vesicle replacement rate + +AN_IHCsynapseParams.x= 10; % replenishment from re-uptake store +AN_IHCsynapseParams.x= 60; % replenishment from re-uptake store + +% reduce l to increase saturated rate +AN_IHCsynapseParams.l= 150; % *loss rate of vesicles from the cleft +% AN_IHCsynapseParams.l= 250; % *loss rate of vesicles from the cleft + +AN_IHCsynapseParams.r= 500; % *reuptake rate from cleft into cell +% AN_IHCsynapseParams.r= 300; % *reuptake rate from cleft into cell + +AN_IHCsynapseParams.refractory_period= 0.00075; +% number of AN fibers at each BF (used only for spike generation) +AN_IHCsynapseParams.numFibers= 100; +AN_IHCsynapseParams.TWdelay=0.004; % ?delay before stimulus first spike + +AN_IHCsynapseParams.ANspeedUpFactor=5; % longer epochs for computing spikes. + +%% #7 MacGregorMulti (first order brainstem neurons) +MacGregorMultiParams=[]; +MacGregorMultiType='chopper'; % MacGregorMultiType='primary-like'; %choose +switch MacGregorMultiType + case 'primary-like' + MacGregorMultiParams.nNeuronsPerBF= 10; % N neurons per BF + MacGregorMultiParams.type = 'primary-like cell'; + MacGregorMultiParams.fibersPerNeuron=4; % N input fibers + MacGregorMultiParams.dendriteLPfreq=200; % dendritic filter + MacGregorMultiParams.currentPerSpike=0.11e-6; % (A) per spike + MacGregorMultiParams.Cap=4.55e-9; % cell capacitance (Siemens) + MacGregorMultiParams.tauM=5e-4; % membrane time constant (s) + MacGregorMultiParams.Ek=-0.01; % K+ eq. potential (V) + MacGregorMultiParams.dGkSpike=3.64e-5; % K+ cond.shift on spike,S + MacGregorMultiParams.tauGk= 0.0012; % K+ conductance tau (s) + MacGregorMultiParams.Th0= 0.01; % equilibrium threshold (V) + MacGregorMultiParams.c= 0.01; % threshold shift on spike, (V) + MacGregorMultiParams.tauTh= 0.015; % variable threshold tau + MacGregorMultiParams.Er=-0.06; % resting potential (V) + MacGregorMultiParams.Eb=0.06; % spike height (V) + + case 'chopper' + MacGregorMultiParams.nNeuronsPerBF= 10; % N neurons per BF + MacGregorMultiParams.type = 'chopper cell'; + MacGregorMultiParams.fibersPerNeuron=10; % N input fibers +% MacGregorMultiParams.fibersPerNeuron=6; % N input fibers + + MacGregorMultiParams.dendriteLPfreq=50; % dendritic filter + MacGregorMultiParams.currentPerSpike=35e-9; % *per spike + MacGregorMultiParams.currentPerSpike=30e-9; % *per spike + + MacGregorMultiParams.Cap=1.67e-8; % ??cell capacitance (Siemens) + MacGregorMultiParams.tauM=0.002; % membrane time constant (s) + MacGregorMultiParams.Ek=-0.01; % K+ eq. potential (V) + MacGregorMultiParams.dGkSpike=1.33e-4; % K+ cond.shift on spike,S + MacGregorMultiParams.tauGk= 0.0005;% K+ conductance tau (s) + MacGregorMultiParams.Th0= 0.01; % equilibrium threshold (V) + MacGregorMultiParams.c= 0; % threshold shift on spike, (V) + MacGregorMultiParams.tauTh= 0.02; % variable threshold tau + MacGregorMultiParams.Er=-0.06; % resting potential (V) + MacGregorMultiParams.Eb=0.06; % spike height (V) + MacGregorMultiParams.PSTHbinWidth= 1e-4; +end + +%% #8 MacGregor (second-order neuron). Only one per channel +MacGregorParams=[]; % clear the structure first +MacGregorParams.type = 'chopper cell'; +MacGregorParams.fibersPerNeuron=10; % N input fibers +MacGregorParams.dendriteLPfreq=100; % dendritic filter +MacGregorParams.currentPerSpike=120e-9;% *(A) per spike +MacGregorParams.currentPerSpike=30e-9;% *(A) per spike + +MacGregorParams.Cap=16.7e-9; % cell capacitance (Siemens) +MacGregorParams.tauM=0.002; % membrane time constant (s) +MacGregorParams.Ek=-0.01; % K+ eq. potential (V) +MacGregorParams.dGkSpike=1.33e-4; % K+ cond.shift on spike,S +MacGregorParams.tauGk= 0.0005; % K+ conductance tau (s) +MacGregorParams.Th0= 0.01; % equilibrium threshold (V) +MacGregorParams.c= 0; % threshold shift on spike, (V) +MacGregorParams.tauTh= 0.02; % variable threshold tau +MacGregorParams.Er=-0.06; % resting potential (V) +MacGregorParams.Eb=0.06; % spike height (V) +MacGregorParams.debugging=0; % (special) +% wideband accepts input from all channels (of same fiber type) +% use wideband to create inhibitory units +MacGregorParams.wideband=0; % special for wideband units +% MacGregorParams.saveAllData=0; + +%% #9 filteredSACF +minPitch= 300; maxPitch= 3000; numPitches=60; % specify lags +pitches=100*log10(logspace(minPitch/100, maxPitch/100, numPitches)); +filteredSACFParams.lags=1./pitches; % autocorrelation lags vector +filteredSACFParams.acfTau= .003; % time constant of running ACF +filteredSACFParams.lambda= 0.12; % slower filter to smooth ACF +filteredSACFParams.plotFilteredSACF=1; % 0 plots unfiltered ACFs +filteredSACFParams.plotACFs=0; % special plot (see code) +% filteredSACFParams.usePressnitzer=0; % attenuates ACF at long lags +filteredSACFParams.lagsProcedure= 'useAllLags'; +% filteredSACFParams.lagsProcedure= 'useBernsteinLagWeights'; +% filteredSACFParams.lagsProcedure= 'omitShortLags'; +filteredSACFParams.criterionForOmittingLags=3; + +% checks +if AN_IHCsynapseParams.numFibers<MacGregorMultiParams.fibersPerNeuron + error('MacGregorMulti: too few input fibers for input to MacG unit') +end + + +%% write all parameters to the command window +% showParams is currently set at the top of htis function +if showParams + fprintf('\n %%%%%%%%\n') + fprintf('\n%s\n', method.parameterSource) + fprintf('\n') + nm=UTIL_paramsList(whos); + for i=1:length(nm) + % eval(['UTIL_showStruct(' nm{i} ', ''' nm{i} ''')']) + if ~strcmp(nm(i), 'method') + eval(['UTIL_showStructureSummary(' nm{i} ', ''' nm{i} ''', 10)']) + end + end +end + + + +% ********************************************************************** comparison data +% store individual data here for display on the multiThreshold GUI (if used) +% the final value in each vector is an identifier (BF or duration)) +if isstruct(experiment) + switch experiment.paradigm + case {'IFMC','IFMC_8ms'} + % based on MPa + comparisonData=[ + 66 51 49 48 46 45 54 250; + 60 54 46 42 39 49 65 500; + 64 51 38 32 33 59 75 1000; + 59 51 36 30 41 81 93 2000; + 71 63 53 44 36 76 95 4000; + 70 64 43 35 35 66 88 6000; + 110 110 110 110 110 110 110 8000; + ]; + if length(BFlist)==1 && ~isempty(comparisonData) + availableFrequencies=comparisonData(:,end)'; + findRow= find(BFlist==availableFrequencies); + if ~isempty (findRow) + experiment.comparisonData=comparisonData(findRow,:); + end + end + + case {'TMC','TMC_8ms'} + % based on MPa + comparisonData=[ + 48 58 63 68 75 80 85 92 99 250; + 33 39 40 49 52 61 64 77 79 500; + 39 42 50 81 83 92 96 97 110 1000; + 24 26 32 37 46 51 59 71 78 2000; + 65 68 77 85 91 93 110 110 110 4000; + 20 19 26 44 80 95 96 110 110 6000; + ]; + if length(BFlist)==1 && ~isempty(comparisonData) + availableFrequencies=comparisonData(:,end)'; + findRow= find(BFlist==availableFrequencies); + if ~isempty (findRow) + experiment.comparisonData=comparisonData(findRow,:); + end + end + + case { 'absThreshold', 'absThreshold_8'} + % MPa thresholds + experiment.comparisonData=[ + 32 26 16 18 22 22 0.008; + 16 13 6 9 15 11 0.500 + ]; + + + otherwise + experiment.comparisonData=[]; + end +end + +
--- a/parameterStore/MAPparamsOHC.m Fri Jun 17 14:30:28 2011 +0100 +++ b/parameterStore/MAPparamsOHC.m Tue Jun 21 14:58:12 2011 +0100 @@ -1,5 +1,5 @@ function method=MAPparamsOHC ... - (BFlist, sampleRate, showParams) + (BFlist, sampleRate, showParams, paramChanges) % MAPparams<> establishes a complete set of MAP parameters % Parameter file names must be of the form <MAPparams> <name> % @@ -58,11 +58,11 @@ % i.e. a minimum ratio of 0.056. % 'spikes' model: AR based on brainstem spiking activity (LSR) OMEParams.rateToAttenuationFactor=0.006; % * N(all ICspikes) -% OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) +% OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) % 'probability model': Ar based on AN firing probabilities (LSR) OMEParams.rateToAttenuationFactorProb=0.01;% * N(all ANrates) -% OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates) +% OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates) % asymptote should be around 100-200 ms OMEParams.ARtau=.05; % AR smoothing function @@ -104,9 +104,9 @@ DRNLParams.rateToAttenuationFactor = .01; % strength of MOC % DRNLParams.rateToAttenuationFactor = 0; % strength of MOC % 'probability' model: MOC based on AN spiking activity (HSR) -DRNLParams.rateToAttenuationFactorProb = .005; % strength of MOC +DRNLParams.rateToAttenuationFactorProb = .0055; % strength of MOC % DRNLParams.rateToAttenuationFactorProb = .0; % strength of MOC -DRNLParams.MOCrateThreshold =70; % spikes/s probability only +DRNLParams.MOCrateThresholdProb =70; % spikes/s probability only DRNLParams.MOCtau =.1; % smoothing for MOC @@ -263,6 +263,16 @@ end +%% now accept last minute parameter changes required by the calling program +% paramChanges +if nargin>3 && ~isempty(paramChanges) + nChanges=length(paramChanges); + for idx=1:nChanges + eval(paramChanges{idx}) + end +end + + %% write all parameters to the command window % showParams is currently set at the top of htis function if showParams @@ -276,6 +286,14 @@ eval(['UTIL_showStructureSummary(' nm{i} ', ''' nm{i} ''', 10)']) end end + + % highlight parameter changes made locally + if nargin>3 && ~isempty(paramChanges) + fprintf('\n Local parameter changes:\n') + for i=1:length(paramChanges) + disp(paramChanges{i}) + end + end end
--- a/parameterStore/MAPparamsStd.m Fri Jun 17 14:30:28 2011 +0100 +++ b/parameterStore/MAPparamsStd.m Tue Jun 21 14:58:12 2011 +0100 @@ -117,7 +117,7 @@ DRNLParams.rateToAttenuationFactorProb = 0; % 0 = MOC off (spikes) end DRNLParams.MOCtau =.03; % smoothing for MOC -DRNLParams.MOCrateThreshold =50; % set to AN rate threshold +DRNLParams.MOCrateThresholdProb =50; % set to AN rate threshold %% #4 IHC_cilia_RPParams
--- a/testPrograms/demoTwisterProbability.m Fri Jun 17 14:30:28 2011 +0100 +++ b/testPrograms/demoTwisterProbability.m Tue Jun 21 14:58:12 2011 +0100 @@ -63,18 +63,23 @@ %% run the model tic +fprintf('\n') +disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)]) +disp([num2str(numChannels) ' channel model']) +disp('Computing ...') + MAP1_14(inputSignal, sampleRate, BFlist, ... MAPparamsName, AN_spikesOrProbability, paramChanges); -toc % the model run is now complete. Now display the results -UTIL_showMAP(showMapOptions) +UTIL_showMAP(showMapOptions, paramChanges) toc path(restorePath) + function inputSignal=createMultiTone(sampleRate, toneFrequency, ... leveldBSPL, duration, rampDuration) % Create pure tone stimulus
--- a/testPrograms/demoTwisterSpikes.m Fri Jun 17 14:30:28 2011 +0100 +++ b/testPrograms/demoTwisterSpikes.m Tue Jun 21 14:58:12 2011 +0100 @@ -45,9 +45,9 @@ showMapOptions.printFiringRates=1; showMapOptions.showACF=0; showMapOptions.showEfferent=0; - + showMapOptions.surfSpikes=1; + %% Generate stimuli - switch signalType case 'tones' inputSignal=createMultiTone(sampleRate, toneFrequency, ... @@ -70,17 +70,17 @@ disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)]) disp([num2str(numChannels) ' channel model']) disp('Computing ...') + MAP1_14(inputSignal, sampleRate, BFlist, ... MAPparamsName, AN_spikesOrProbability, paramChanges); -toc + % the model run is now complete. Now display the results -UTIL_showMAP(showMapOptions) +UTIL_showMAP(showMapOptions, paramChanges) toc path(restorePath) - function inputSignal=createMultiTone(sampleRate, toneFrequency, ... leveldBSPL, duration, rampDuration) % Create pure tone stimulus
--- a/testPrograms/test_MAP1_14.m Fri Jun 17 14:30:28 2011 +0100 +++ b/testPrograms/test_MAP1_14.m Tue Jun 21 14:58:12 2011 +0100 @@ -42,29 +42,30 @@ %% #2 probability (fast) or spikes (slow) representation -AN_spikesOrProbability='spikes'; +% AN_spikesOrProbability='spikes'; % or % NB probabilities are not corrected for refractory effects -% AN_spikesOrProbability='probability'; +AN_spikesOrProbability='probability'; %% #3 pure tone, harmonic sequence or speech file input signalType= 'tones'; -sampleRate= 100000; -duration=0.050; % seconds +sampleRate= 50000; +duration=0.250; % seconds % toneFrequency= 250:250:8000; % harmonic sequence (Hz) toneFrequency= 1000; % or a pure tone (Hz8 rampDuration=.005; % seconds % or + % signalType= 'file'; % fileName='twister_44kHz'; %% #4 rms level % signal details -leveldBSPL= 70; % dB SPL +leveldBSPL= 30; % dB SPL %% #5 number of channels in the model @@ -90,10 +91,10 @@ % paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ... % 'IHCpreSynapseParams.tauCa=86e-6;'}; -% paramChanges={'DRNLParams.rateToAttenuationFactorProb = 0;'}; +paramChanges={'DRNLParams.rateToAttenuationFactorProb = 0;'}; % paramChanges={'IHCpreSynapseParams.tauCa=86e-6;', -% 'AN_IHCsynapseParams.numFibers= 100;'}; +% 'AN_IHCsynapseParams.numFibers= 1000;'}; %% delare 'showMap' options to control graphical output @@ -152,15 +153,10 @@ MAP1_14(inputSignal, sampleRate, BFlist, ... MAPparamsName, AN_spikesOrProbability, paramChanges); -toc + % the model run is now complete. Now display the results -disp(' param changes to list of parameters below') -for i=1:length(paramChanges) - disp(paramChanges{i}) -end -UTIL_showMAP(showMapOptions) - +UTIL_showMAP(showMapOptions, paramChanges) toc path(restorePath) @@ -187,6 +183,5 @@ % add 10 ms silence silence= zeros(1,round(0.05/dt)); -silence= zeros(1,round(0.05/dt)); inputSignal= [silence inputSignal silence];
--- a/utilities/UTIL_PSTHmaker.m Fri Jun 17 14:30:28 2011 +0100 +++ b/utilities/UTIL_PSTHmaker.m Tue Jun 21 14:58:12 2011 +0100 @@ -8,9 +8,7 @@ % inputData is a channel x time matrix % PSTH is the reduced matrix, the sum of all elements in the bin % -% mandatory fileds: -% method.dt -% method.PSTHbinWidth + [numChannels numDataPoints]= size(inputData); % Multiple fibers is the same as repeat trials
--- a/utilities/UTIL_showMAP.m Fri Jun 17 14:30:28 2011 +0100 +++ b/utilities/UTIL_showMAP.m Tue Jun 21 14:58:12 2011 +0100 @@ -1,4 +1,4 @@ -function UTIL_showMAP (showMapOptions) +function UTIL_showMAP (showMapOptions, paramChanges) % UTIL_showMAP produces summaries of the output from MAP's mmost recent run % All MAP outputs are stored in global variables and UTIL_showMAP % simply assumes that they are in place. @@ -40,15 +40,21 @@ if ~isfield(showMapOptions,'fileName'),showMapOptions.fileName=[];end if ~isfield(showMapOptions,'surfSpikes'),showMapOptions.surfSpikes=0;end - +%% send all model parameters to command window if showMapOptions.printModelParameters % Read parameters from MAPparams<***> file in 'parameterStore' folder % and print out all parameters - cmd=['MAPparams' saveMAPparamsName ... - '(-1, 1/dt, 1);']; - eval(cmd); + if nargin>1 + cmd=['MAPparams' saveMAPparamsName '(-1, 1/dt, 1, paramChanges);']; + eval(cmd); + else + cmd=['MAPparams' saveMAPparamsName '(-1, 1/dt, 1);']; + eval(cmd); + disp(' no parameter changes found') + end end +%% summarise firing rates in command window if showMapOptions.printFiringRates %% print summary firing rates fprintf('\n\n') @@ -71,14 +77,14 @@ % disp(['IC by type: ' num2str(mean(ICfiberTypeRates,2)')]) else disp(['AN: ' num2str(mean(mean(ANprobRateOutput)))]) - [PSTH pointsPerBin]= UTIL_makePSTH(ANprobRateOutput, dt, 0.001); + PSTH= UTIL_PSTHmakerb(ANprobRateOutput, dt, 0.001); disp(['max max AN: ' num2str(max(max(... - PSTH/pointsPerBin )))]) + PSTH )))]) end end -%% figure (99) summarises main model output +%% figure (99): display output from all model stages if showMapOptions.showModelOutput plotInstructions.figureNo=99; signalRMS=mean(savedInputSignal.^2)^0.5; @@ -150,8 +156,10 @@ UTIL_plotMatrix(ICmembraneOutput, plotInstructions); end - otherwise % probability (4-6) - plotInstructions.displaydt=dt; + otherwise % AN rate based on probability of firing + PSTHbinWidth=0.001; + PSTH= UTIL_PSTHmakerb(ANprobRateOutput, dt, PSTHbinWidth); + plotInstructions.displaydt=PSTHbinWidth; plotInstructions.numPlots=2; plotInstructions.subPlotNo=2; plotInstructions.yLabel='BF'; @@ -159,13 +167,15 @@ plotInstructions.yLabel='LSR HSR'; plotInstructions.plotDivider=1; end - plotInstructions.title='AN - spike probability'; - UTIL_plotMatrix(ANprobRateOutput, plotInstructions); + plotInstructions.title='AN - spike rate'; + UTIL_plotMatrix(PSTH, plotInstructions); end end if showMapOptions.surfProbability %% surface plot of probability + if strcmp(saveAN_spikesOrProbability,'probability') && ... + length(savedBFlist)>2 figure(97), clf % select only HSR fibers at the bottom of the matrix ANprobRateOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:); @@ -185,15 +195,16 @@ view([-20 60]) % view([0 90]) title ([showMapOptions.fileName ': ' num2str(signalRMSdb,'% 3.0f') ' dB']) + end end if showMapOptions.surfSpikes - %% surface plot of probability + %% surface plot of AN spikes figure(97), clf % select only HSR fibers at the bottom of the matrix ANoutput= ANoutput(end-length(savedBFlist)+1:end,:); - PSTHbinWidth=0.001; % 1 ms bins - PSTH=UTIL_PSTHmakerb(ANoutput, ANdt, PSTHbinWidth); + PSTHbinWidth=0.005; % 1 ms bins + PSTH=UTIL_makePSTH(ANoutput, ANdt, PSTHbinWidth); [nY nX]=size(PSTH); time=PSTHbinWidth*(1:nX); surf(time, savedBFlist, PSTH) @@ -211,7 +222,7 @@ end -%% plot efferent control values as dB +%% figure(98) plot efferent control values as dB if showMapOptions.showEfferent plotInstructions=[]; plotInstructions.figureNo=98;