Mercurial > hg > map
changeset 16:37a379b27cff
changes to showMap
| author | Ray Meddis <rmeddis@essex.ac.uk> |
|---|---|
| date | Tue, 07 Jun 2011 09:53:50 +0100 |
| parents | 35af36fe0a35 |
| children | e9e263e4fcde c489ebada16e |
| files | MAP/temp.m multithreshold 1.46/savedData/mostRecentResults.mat parameterStore/MAPparamsEndo.m parameterStore/MAPparamsNormal.m parameterStore/MAPparamsOHC.m testPrograms/showMAP.m testPrograms/test_MAP1_14.m wavFileStore/babble.wav |
| diffstat | 8 files changed, 1247 insertions(+), 85 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MAP/temp.m Tue Jun 07 09:53:50 2011 +0100 @@ -0,0 +1,1128 @@ + +function MAP1_14(inputSignal, sampleRate, BFlist, MAPparamsName, ... + AN_spikesOrProbability, paramChanges) +% To test this function use test_MAP1_14 in this folder +% +% All arguments are mandatory. +% +% BFlist is a list of BFs but can be '-1' to allow MAPparams to choose +% + +% MAPparamsName='Normal'; % source of model parameters +% AN_spikesOrProbability='spikes'; % or 'probability' +% paramChanges is a cell array of strings that can be used to make last +% minute parameter changes, e.g., to simulate OHC loss +% paramChanges{1}= 'DRNLParams.a=0;'; + +% The model parameters are established in the MAPparams<***> file +% and stored as global + +restorePath=path; +addpath (['..' filesep 'parameterStore']) + +global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams +global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams + +% All of the results of this function are stored as global +global dt ANdt savedBFlist saveAN_spikesOrProbability saveMAPparamsName... + savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ... + DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV... + IHCoutput ANprobRateOutput ANoutput savePavailable tauCas ... + CNoutput ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation + +% Normally only ICoutput(logical spike matrix) or ANprobRateOutput will be +% needed by the user; so the following will suffice +% global ANdt ICoutput ANprobRateOutput + +% Note that sampleRate has not changed from the original function call and +% ANprobRateOutput is sampled at this rate +% However ANoutput, CNoutput and IC output are stored as logical +% 'spike' matrices using a lower sample rate (see ANdt). + +% When AN_spikesOrProbability is set to probability, +% no spike matrices are computed. +% When AN_spikesOrProbability is set to 'spikes', +% no probability output is computed + +% Efferent control variables are ARattenuation and MOCattenuation +% These are scalars between 1 (no attenuation) and 0. +% They are represented with dt=1/sampleRate (not ANdt) +% They are computed using either AN probability rate output +% or IC (spikes) output as approrpriate. +% AR is computed using across channel activity +% MOC is computed on a within-channel basis. + + +% 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) +segmentDuration=method.segmentDuration; +segmentLength=round(segmentDuration/ dt); +segmentTime=dt*(1:segmentLength); % used in debugging plots + +% all spiking activity is computed using longer epochs +ANspeedUpFactor=AN_IHCsynapseParams.ANspeedUpFactor; % e.g.5 times + +% inputSignal must be row vector +[r c]=size(inputSignal); +if r>c, inputSignal=inputSignal'; end % transpose +% ignore stereo signals +inputSignal=inputSignal(1,:); % drop any second channel +savedInputSignal=inputSignal; + +% Segment the signal +% The sgment length is given but the signal length must be adjusted to be a +% multiple of both the segment length and the reduced segmentlength +[nSignalRows signalLength]=size(inputSignal); +segmentLength=ceil(segmentLength/ANspeedUpFactor)*ANspeedUpFactor; +% Make the signal length a whole multiple of the segment length +nSignalSegments=ceil(signalLength/segmentLength); +padSize=nSignalSegments*segmentLength-signalLength; +pad=zeros(nSignalRows,padSize); +inputSignal=[inputSignal pad]; +[ignore signalLength]=size(inputSignal); + +% AN (spikes) is computed at a lower sample rate when spikes required +% so it has a reduced segment length (see 'ANspeeUpFactor' above) +% AN CN and IC all use this sample interval +ANdt=dt*ANspeedUpFactor; +reducedSegmentLength=round(segmentLength/ANspeedUpFactor); +reducedSignalLength= round(signalLength/ANspeedUpFactor); + +%% Initialise with respect to each stage before computing +% by allocating memory, +% by computing constants +% by establishing easy to read variable names +% The computations are made in segments and boundary conditions must +% be established and stored. These are found in variables with +% 'boundary' or 'bndry' in the name + +%% OME --- +% external ear resonances +OMEexternalResonanceFilters=OMEParams.externalResonanceFilters; +[nOMEExtFilters c]=size(OMEexternalResonanceFilters); +% details of external (outer ear) resonances +OMEgaindBs=OMEexternalResonanceFilters(:,1); +OMEgainScalars=10.^(OMEgaindBs/20); +OMEfilterOrder=OMEexternalResonanceFilters(:,2); +OMElowerCutOff=OMEexternalResonanceFilters(:,3); +OMEupperCutOff=OMEexternalResonanceFilters(:,4); +% external resonance coefficients +ExtFilter_b=cell(nOMEExtFilters,1); +ExtFilter_a=cell(nOMEExtFilters,1); +for idx=1:nOMEExtFilters + Nyquist=sampleRate/2; + [b, a] = butter(OMEfilterOrder(idx), ... + [OMElowerCutOff(idx) OMEupperCutOff(idx)]... + /Nyquist); + ExtFilter_b{idx}=b; + ExtFilter_a{idx}=a; +end +OMEExtFilterBndry=cell(2,1); +OMEextEarPressure=zeros(1,signalLength); % pressure at tympanic membrane + +% pressure to velocity conversion using smoothing filter (50 Hz cutoff) +tau=1/(2*pi*50); +a1=dt/tau-1; a0=1; +b0=1+ a1; +TMdisp_b=b0; TMdisp_a=[a0 a1]; +% figure(9), freqz(TMdisp_b, TMdisp_a) +OME_TMdisplacementBndry=[]; + +% OME high pass (simulates poor low frequency stapes response) +OMEhighPassHighCutOff=OMEParams.OMEstapesLPcutoff; +Nyquist=sampleRate/2; +[stapesDisp_b,stapesDisp_a] = butter(1, OMEhighPassHighCutOff/Nyquist, 'high'); +% figure(10), freqz(stapesDisp_b, stapesDisp_a) + +OMEhighPassBndry=[]; + +% OMEampStapes might be reducdant (use OMEParams.stapesScalar) +stapesScalar= OMEParams.stapesScalar; + +% Acoustic reflex +efferentDelayPts=round(OMEParams.ARdelay/dt); +% smoothing filter +a1=dt/OMEParams.ARtau-1; a0=1; +b0=1+ a1; +ARfilt_b=b0; ARfilt_a=[a0 a1]; + +ARattenuation=ones(1,signalLength); +ARrateThreshold=OMEParams.ARrateThreshold; % may not be used +ARrateToAttenuationFactor=OMEParams.rateToAttenuationFactor; +ARrateToAttenuationFactorProb=OMEParams.rateToAttenuationFactorProb; +ARboundary=[]; +ARboundaryProb=0; + +% save complete OME record (stapes displacement) +OMEoutput=zeros(1,signalLength); +TMoutput=zeros(1,signalLength); + +%% BM --- +% BM is represented as a list of locations identified by BF +DRNL_BFs=BFlist; +nBFs= length(DRNL_BFs); + +% DRNLchannelParameters=DRNLParams.channelParameters; +DRNLresponse= zeros(nBFs, segmentLength); + +MOCrateToAttenuationFactor=DRNLParams.rateToAttenuationFactor; +rateToAttenuationFactorProb=DRNLParams.rateToAttenuationFactorProb; +MOCrateThreshold=DRNLParams.MOCrateThreshold; + +% smoothing filter for MOC +a1=dt/DRNLParams.MOCtau-1; a0=1; +b0=1+ a1; +MOCfilt_b=b0; MOCfilt_a=[a0 a1]; +% figure(9), freqz(stapesDisp_b, stapesDisp_a) +MOCboundary=cell(nBFs,1); +MOCprobBoundary=cell(nBFs,1); + +MOCattSegment=zeros(nBFs,reducedSegmentLength); +MOCattenuation=ones(nBFs,signalLength); + +if DRNLParams.a>0 + DRNLcompressionThreshold=10^((1/(1-DRNLParams.c))* ... + log10(DRNLParams.b/DRNLParams.a)); +else + DRNLcompressionThreshold=inf; +end + +DRNLlinearOrder= DRNLParams.linOrder; +DRNLnonlinearOrder= DRNLParams.nonlinOrder; + +DRNLa=DRNLParams.a; +DRNLb=DRNLParams.b; +DRNLc=DRNLParams.c; +linGAIN=DRNLParams.g; +% +% gammatone filter coefficients for linear pathway +bw=DRNLParams.linBWs'; +phi = 2 * pi * bw * dt; +cf=DRNLParams.linCFs'; +theta = 2 * pi * cf * dt; +cos_theta = cos(theta); +sin_theta = sin(theta); +alpha = -exp(-phi).* cos_theta; +b0 = ones(nBFs,1); +b1 = 2 * alpha; +b2 = exp(-2 * phi); +z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i; +z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i; +z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i; +tf = (z2 + z3) ./ z1; +a0 = abs(tf); +a1 = alpha .* a0; +GTlin_a = [b0, b1, b2]; +GTlin_b = [a0, a1]; +GTlinOrder=DRNLlinearOrder; +GTlinBdry=cell(nBFs,GTlinOrder); + +% nonlinear gammatone filter coefficients +bw=DRNLParams.nlBWs'; +phi = 2 * pi * bw * dt; +cf=DRNLParams.nonlinCFs'; +theta = 2 * pi * cf * dt; +cos_theta = cos(theta); +sin_theta = sin(theta); +alpha = -exp(-phi).* cos_theta; +b0 = ones(nBFs,1); +b1 = 2 * alpha; +b2 = exp(-2 * phi); +z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i; +z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i; +z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i; +tf = (z2 + z3) ./ z1; +a0 = abs(tf); +a1 = alpha .* a0; +GTnonlin_a = [b0, b1, b2]; +GTnonlin_b = [a0, a1]; +GTnonlinOrder=DRNLnonlinearOrder; +GTnonlinBdry1=cell(nBFs, GTnonlinOrder); +GTnonlinBdry2=cell(nBFs, GTnonlinOrder); + +% complete BM record (BM displacement) +DRNLoutput=zeros(nBFs, signalLength); + + +%% IHC --- +% IHC cilia activity and receptor potential +% viscous coupling between BM and stereocilia displacement +% Nyquist=sampleRate/2; +% IHCcutoff=1/(2*pi*IHC_cilia_RPParams.tc); +% [IHCciliaFilter_b,IHCciliaFilter_a]=... +% butter(1, IHCcutoff/Nyquist, 'high'); +a1=dt/IHC_cilia_RPParams.tc-1; a0=1; +b0=1+ a1; +% high pass (i.e. low pass reversed) +IHCciliaFilter_b=[a0 a1]; IHCciliaFilter_a=b0; +% figure(9), freqz(IHCciliaFilter_b, IHCciliaFilter_a) + +IHCciliaBndry=cell(nBFs,1); + +% IHC apical conductance (Boltzman function) +IHC_C= IHC_cilia_RPParams.C; +IHCu0= IHC_cilia_RPParams.u0; +IHCu1= IHC_cilia_RPParams.u1; +IHCs0= IHC_cilia_RPParams.s0; +IHCs1= IHC_cilia_RPParams.s1; +IHCGmax= IHC_cilia_RPParams.Gmax; +IHCGa= IHC_cilia_RPParams.Ga; % (leakage) + +IHCGu0 = IHCGa+IHCGmax./(1+exp(IHCu0/IHCs0).*(1+exp(IHCu1/IHCs1))); +IHCrestingCiliaCond=IHCGu0; + +% Receptor potential +IHC_Cab= IHC_cilia_RPParams.Cab; +IHC_Gk= IHC_cilia_RPParams.Gk; +IHC_Et= IHC_cilia_RPParams.Et; +IHC_Ek= IHC_cilia_RPParams.Ek; +IHC_Ekp= IHC_Ek+IHC_Et*IHC_cilia_RPParams.Rpc; + +IHCrestingV= (IHC_Gk*IHC_Ekp+IHCGu0*IHC_Et)/(IHCGu0+IHC_Gk); + +IHC_Vnow= IHCrestingV*ones(nBFs,1); % initial voltage +IHC_RP= zeros(nBFs,segmentLength); + +% complete record of IHC receptor potential (V) +IHCciliaDisplacement= zeros(nBFs,segmentLength); +IHCoutput= zeros(nBFs,signalLength); +IHC_cilia_output= zeros(nBFs,signalLength); + + +%% pre-synapse --- +% Each BF is replicated using a different fiber type to make a 'channel' +% The number of channels is nBFs x nANfiberTypes +% Fiber types are specified in terms of tauCa +nANfiberTypes= length(IHCpreSynapseParams.tauCa); +tauCas= IHCpreSynapseParams.tauCa; +nChannels= nANfiberTypes*nBFs; +synapticCa= zeros(nChannels,segmentLength); + +% Calcium control (more calcium, greater release rate) +ECa=IHCpreSynapseParams.ECa; +gamma=IHCpreSynapseParams.gamma; +beta=IHCpreSynapseParams.beta; +tauM=IHCpreSynapseParams.tauM; +mICa=zeros(nChannels,segmentLength); +GmaxCa=IHCpreSynapseParams.GmaxCa; +synapse_z= IHCpreSynapseParams.z; +synapse_power=IHCpreSynapseParams.power; + +% tauCa vector is established across channels to allow vectorization +% (one tauCa per channel). Do not confuse with tauCas (one pre fiber type) +tauCa=repmat(tauCas, nBFs,1); +tauCa=reshape(tauCa, nChannels, 1); + +% presynapse startup values (vectors, length:nChannels) +% proportion (0 - 1) of Ca channels open at IHCrestingV +mICaCurrent=((1+beta^-1 * exp(-gamma*IHCrestingV))^-1)... + *ones(nBFs*nANfiberTypes,1); +% corresponding startup currents +ICaCurrent= (GmaxCa*mICaCurrent.^3) * (IHCrestingV-ECa); +CaCurrent= ICaCurrent.*tauCa; + +% vesicle release rate at startup (one per channel) +% kt0 is used only at initialisation +kt0= -synapse_z * CaCurrent.^synapse_power; + + +%% AN --- +% each row of the AN matrices represents one AN fiber +% The results computed either for probabiities *or* for spikes (not both) +% Spikes are necessary if CN and IC are to be computed +nFibersPerChannel= AN_IHCsynapseParams.numFibers; +nANfibers= nChannels*nFibersPerChannel; + +y=AN_IHCsynapseParams.y; +l=AN_IHCsynapseParams.l; +x=AN_IHCsynapseParams.x; +r=AN_IHCsynapseParams.r; +M=round(AN_IHCsynapseParams.M); + +% probability (NB initial 'P' on everything) +PAN_ydt = repmat(AN_IHCsynapseParams.y*dt, nChannels,1); +PAN_ldt = repmat(AN_IHCsynapseParams.l*dt, nChannels,1); +PAN_xdt = repmat(AN_IHCsynapseParams.x*dt, nChannels,1); +PAN_rdt = repmat(AN_IHCsynapseParams.r*dt, nChannels,1); +PAN_rdt_plus_ldt = PAN_rdt + PAN_ldt; +PAN_M=round(AN_IHCsynapseParams.M); + +% compute starting values +Pcleft = kt0* y* M ./ (y*(l+r)+ kt0* l); +Pavailable = Pcleft*(l+r)./kt0; +Preprocess = Pcleft*r/x; % canbe fractional + +ANprobability=zeros(nChannels,segmentLength); +ANprobRateOutput=zeros(nChannels,signalLength); +% 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); +AN_ldt= repmat(AN_IHCsynapseParams.l*ANdt, nANfibers,1); +AN_xdt= repmat(AN_IHCsynapseParams.x*ANdt, nANfibers,1); +AN_rdt= repmat(AN_IHCsynapseParams.r*ANdt, nANfibers,1); +AN_rdt_plus_ldt= AN_rdt + AN_ldt; +AN_M= round(AN_IHCsynapseParams.M); + +% kt0 is initial release rate +% Establish as a vector (length=channel x number of fibers) +kt0= repmat(kt0', nFibersPerChannel, 1); +kt0=reshape(kt0, nANfibers,1); + +% starting values for reservoirs +AN_cleft = kt0* y* M ./ (y*(l+r)+ kt0* l); +AN_available = round(AN_cleft*(l+r)./kt0); %must be integer +AN_reprocess = AN_cleft*r/x; + +% output is in a logical array spikes = 1/ 0. +ANspikes= false(nANfibers,reducedSegmentLength); +ANoutput= false(nANfibers,reducedSignalLength); + + +%% CN (first brain stem nucleus - could be any subdivision of CN) +% Input to a CN neuorn is a random selection of AN fibers within a channel +% The number of AN fibers used is ANfibersFanInToCN +ANfibersFanInToCN=MacGregorMultiParams.fibersPerNeuron; +nCNneuronsPerChannel=MacGregorMultiParams.nNeuronsPerBF; +% CNtauGk (Potassium time constant) determines the rate of firing of +% the unit when driven hard by a DC input (not normally >350 sp/s) +CNtauGk=MacGregorMultiParams.tauGk; +ANavailableFibersPerChan=AN_IHCsynapseParams.numFibers; +nCNneurons=nCNneuronsPerChannel*nChannels; +% nCNneuronsPerFiberType= nCNneurons/nANfiberTypes; + +CNmembranePotential=zeros(nCNneurons,reducedSegmentLength); + +% establish which ANfibers (by name) feed into which CN nuerons +CNinputfiberLists=zeros(nChannels*nCNneuronsPerChannel, ANfibersFanInToCN); +unitNo=1; +for ch=1:nChannels + % Each channel contains a number of units =length(listOfFanInValues) + for idx=1:nCNneuronsPerChannel + fibersUsed=(ch-1)*ANavailableFibersPerChan + ... + ceil(rand(1,ANfibersFanInToCN)* ANavailableFibersPerChan); + CNinputfiberLists(unitNo,:)=fibersUsed; + unitNo=unitNo+1; + end +end + +% input to CN units +AN_PSTH=zeros(nCNneurons,reducedSegmentLength); + +% Generate CNalphaFunction function +% by which spikes are converted to post-synaptic currents +CNdendriteLPfreq= MacGregorMultiParams.dendriteLPfreq; +CNcurrentPerSpike=MacGregorMultiParams.currentPerSpike; +CNspikeToCurrentTau=1/(2*pi*CNdendriteLPfreq); +t=ANdt:ANdt:5*CNspikeToCurrentTau; +CNalphaFunction= (1 / ... + CNspikeToCurrentTau)*t.*exp(-t /CNspikeToCurrentTau); +CNalphaFunction=CNalphaFunction*CNcurrentPerSpike; + +% figure(98), plot(t,CNalphaFunction) +% working memory for implementing convolution + +CNcurrentTemp=... + zeros(nCNneurons,reducedSegmentLength+length(CNalphaFunction)-1); +% trailing alphas are parts of humps carried forward to the next segment +CNtrailingAlphas=zeros(nCNneurons,length(CNalphaFunction)); + +CN_tauM=MacGregorMultiParams.tauM; +CN_tauTh=MacGregorMultiParams.tauTh; +CN_cap=MacGregorMultiParams.Cap; +CN_c=MacGregorMultiParams.c; +CN_b=MacGregorMultiParams.dGkSpike; +CN_Ek=MacGregorMultiParams.Ek; +CN_Eb= MacGregorMultiParams.Eb; +CN_Er=MacGregorMultiParams.Er; +CN_Th0= MacGregorMultiParams.Th0; +CN_E= zeros(nCNneurons,1); +CN_Gk= zeros(nCNneurons,1); +CN_Th= MacGregorMultiParams.Th0*ones(nCNneurons,1); +CN_Eb=CN_Eb.*ones(nCNneurons,1); +CN_Er=CN_Er.*ones(nCNneurons,1); +CNtimeSinceLastSpike=zeros(nCNneurons,1); +% tauGk is the main distinction between neurons +% in fact they are all the same in the standard model +tauGk=repmat(CNtauGk,nChannels*nCNneuronsPerChannel,1); + +CN_PSTH=zeros(nChannels,reducedSegmentLength); +CNoutput=false(nCNneurons,reducedSignalLength); + + +%% MacGregor (IC - second nucleus) -------- +nICcells=nChannels; % one cell per channel + +ICspikeWidth=0.00015; % this may need revisiting +epochsPerSpike=round(ICspikeWidth/ANdt); +if epochsPerSpike<1 + error(['MacGregorMulti: sample rate too low to support ' ... + num2str(ICspikeWidth*1e6) ' microsec spikes']); +end + +% short names +IC_tauM=MacGregorParams.tauM; +IC_tauGk=MacGregorParams.tauGk; +IC_tauTh=MacGregorParams.tauTh; +IC_cap=MacGregorParams.Cap; +IC_c=MacGregorParams.c; +IC_b=MacGregorParams.dGkSpike; +IC_Th0=MacGregorParams.Th0; +IC_Ek=MacGregorParams.Ek; +IC_Eb= MacGregorParams.Eb; +IC_Er=MacGregorParams.Er; + +IC_E=zeros(nICcells,1); +IC_Gk=zeros(nICcells,1); +IC_Th=IC_Th0*ones(nICcells,1); + +% Dendritic filtering, all spikes are replaced by CNalphaFunction functions +ICdendriteLPfreq= MacGregorParams.dendriteLPfreq; +ICcurrentPerSpike=MacGregorParams.currentPerSpike; +ICspikeToCurrentTau=1/(2*pi*ICdendriteLPfreq); +t=ANdt:ANdt:3*ICspikeToCurrentTau; +IC_CNalphaFunction= (ICcurrentPerSpike / ... + ICspikeToCurrentTau)*t.*exp(-t / ICspikeToCurrentTau); +% figure(98), plot(t,IC_CNalphaFunction) + +% working space for implementing alpha function +ICcurrentTemp=... + zeros(nICcells,reducedSegmentLength+length(IC_CNalphaFunction)-1); +ICtrailingAlphas=zeros(nICcells, length(IC_CNalphaFunction)); + +ICfiberTypeRates=zeros(nANfiberTypes,reducedSignalLength); +ICoutput=false(nChannels,reducedSignalLength); + +ICmembranePotential=zeros(nICcells,reducedSegmentLength); +ICmembraneOutput=zeros(nICcells,signalLength); + + +%% Main program %% %% %% %% %% %% %% %% %% %% %% %% %% %% + +% Compute the entire model for each segment +segmentStartPTR=1; +reducedSegmentPTR=1; % when sampling rate is reduced +while segmentStartPTR<signalLength + segmentEndPTR=segmentStartPTR+segmentLength-1; + % shorter segments after speed up. + shorterSegmentEndPTR=reducedSegmentPTR+reducedSegmentLength-1; + + iputPressureSegment=inputSignal... + (:,segmentStartPTR:segmentStartPTR+segmentLength-1); + + % segment debugging plots + % figure(98) + % plot(segmentTime,iputPressureSegment), title('signalSegment') + + + % OME ---------------------- + + % OME Stage 1: external resonances. Add to inputSignal pressure wave + y=iputPressureSegment; + for n=1:nOMEExtFilters + % any number of resonances can be used + [x OMEExtFilterBndry{n}] = ... + filter(ExtFilter_b{n},ExtFilter_a{n},... + iputPressureSegment, OMEExtFilterBndry{n}); + x= x* OMEgainScalars(n); + % This is a parallel resonance so add it + y=y+x; + end + iputPressureSegment=y; + OMEextEarPressure(segmentStartPTR:segmentEndPTR)= iputPressureSegment; + + % OME stage 2: convert input pressure (velocity) to + % tympanic membrane(TM) displacement using low pass filter + [TMdisplacementSegment OME_TMdisplacementBndry] = ... + filter(TMdisp_b,TMdisp_a,iputPressureSegment, ... + OME_TMdisplacementBndry); + % and save it + TMoutput(segmentStartPTR:segmentEndPTR)= TMdisplacementSegment; + + % OME stage 3: middle ear high pass effect to simulate stapes inertia + [stapesDisplacement OMEhighPassBndry] = ... + filter(stapesDisp_b,stapesDisp_a,TMdisplacementSegment, ... + OMEhighPassBndry); + + % OME stage 4: apply stapes scalar + stapesDisplacement=stapesDisplacement*stapesScalar; + + % OME stage 5: acoustic reflex stapes attenuation + % Attenuate the TM response using feedback from LSR fiber activity + if segmentStartPTR>efferentDelayPts + stapesDisplacement= stapesDisplacement.*... + ARattenuation(segmentStartPTR-efferentDelayPts:... + segmentEndPTR-efferentDelayPts); + end + + % segment debugging plots + % figure(98) + % plot(segmentTime, stapesDisplacement), title ('stapesDisplacement') + + % and save + OMEoutput(segmentStartPTR:segmentEndPTR)= stapesDisplacement; + + + %% BM ------------------------------ + % Each location is computed separately + for BFno=1:nBFs + + % *linear* path + linOutput = stapesDisplacement * linGAIN; % linear gain + for order = 1 : GTlinOrder + [linOutput GTlinBdry{BFno,order}] = ... + filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, GTlinBdry{BFno,order}); + end + + % *nonLinear* path + % efferent attenuation (0 <> 1) + if segmentStartPTR>efferentDelayPts + MOC=MOCattenuation(BFno, segmentStartPTR-efferentDelayPts:... + segmentEndPTR-efferentDelayPts); + else % no MOC available yet + MOC=ones(1, segmentLength); + end + + % first gammatone filter + for order = 1 : GTnonlinOrder + [nonlinOutput GTnonlinBdry1{BFno,order}] = ... + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ... + stapesDisplacement, GTnonlinBdry1{BFno,order}); + end + + % broken stick instantaneous compression + % nonlinear gain is weakend by MOC before applied to BM response + y= nonlinOutput.*(MOC* DRNLa); % linear section. + % compress those parts of the signal above the compression + % threshold + abs_x = abs(nonlinOutput); + idx=find(abs_x>DRNLcompressionThreshold); + if ~isempty(idx)>0 + y(idx)=sign(nonlinOutput(idx)).*... + (DRNLb*abs_x(idx).^DRNLc); + end + nonlinOutput=y; + + % second filter removes distortion products + for order = 1 : GTnonlinOrder + [ nonlinOutput GTnonlinBdry2{BFno,order}] = ... + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), nonlinOutput, GTnonlinBdry2{BFno,order}); + end + + % combine the two paths to give the DRNL displacement + DRNLresponse(BFno,:)=linOutput+nonlinOutput; + end % BF + + % segment debugging plots + % figure(98) + % if size(DRNLresponse,1)>3 + % imagesc(DRNLresponse) % matrix display + % title('DRNLresponse'); % single or double channel response + % else + % plot(segmentTime, DRNLresponse) + % end + + % and save it + DRNLoutput(:, segmentStartPTR:segmentEndPTR)= DRNLresponse; + + + %% IHC ------------------------------------ + % BM displacement to IHCciliaDisplacement is a high-pass filter + % because of viscous coupling + for idx=1:nBFs + [IHCciliaDisplacement(idx,:) IHCciliaBndry{idx}] = ... + filter(IHCciliaFilter_b,IHCciliaFilter_a, ... + DRNLresponse(idx,:), IHCciliaBndry{idx}); + end + + % apply scalar + IHCciliaDisplacement=IHCciliaDisplacement* IHC_C; + + % and save it + IHC_cilia_output(:,segmentStartPTR:segmentStartPTR+segmentLength-1)=... + IHCciliaDisplacement; + + % compute apical conductance + G=IHCGmax./(1+exp(-(IHCciliaDisplacement-IHCu0)/IHCs0).*... + (1+exp(-(IHCciliaDisplacement-IHCu1)/IHCs1))); + Gu=G + IHCGa; + + % Compute receptor potential + for idx=1:segmentLength + IHC_Vnow=IHC_Vnow+ (-Gu(:, idx).*(IHC_Vnow-IHC_Et)-... + IHC_Gk*(IHC_Vnow-IHC_Ekp))* dt/IHC_Cab; + IHC_RP(:,idx)=IHC_Vnow; + end + + % segment debugging plots + % if size(IHC_RP,1)>3 + % surf(IHC_RP), shading interp, title('IHC_RP') + % else + % plot(segmentTime, IHC_RP) + % end + + % and save it + IHCoutput(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=IHC_RP; + + + %% synapse ----------------------------- + % Compute the vesicle release rate for each fiber type at each BF + % replicate IHC_RP for each fiber type + Vsynapse=repmat(IHC_RP, nANfiberTypes,1); + + % look-up table of target fraction channels open for a given IHC_RP + mICaINF= 1./( 1 + exp(-gamma * Vsynapse) /beta); + % fraction of channel open - apply time constant + for idx=1:segmentLength + % mICaINF is the current 'target' value of mICa + mICaCurrent=mICaCurrent+(mICaINF(:,idx)-mICaCurrent)*dt./tauM; + mICa(:,idx)=mICaCurrent; + end + + ICa= (GmaxCa* mICa.^3) .* (Vsynapse- ECa); + + for idx=1:segmentLength + CaCurrent=CaCurrent + ICa(:,idx)*dt - CaCurrent*dt./tauCa; + synapticCa(:,idx)=CaCurrent; + end + synapticCa=-synapticCa; % treat IHCpreSynapseParams as positive substance + + % NB vesicleReleaseRate is /s and is independent of dt + vesicleReleaseRate = synapse_z * synapticCa.^synapse_power; % rate + + % segment debugging plots + % if size(vesicleReleaseRate,1)>3 + % surf(vesicleReleaseRate), shading interp, title('vesicleReleaseRate') + % else + % plot(segmentTime, vesicleReleaseRate) + % end + + + %% AN + switch AN_spikesOrProbability + case 'probability' + % No refractory effect is applied + for t = 1:segmentLength; + M_Pq=PAN_M-Pavailable; + M_Pq(M_Pq<0)=0; + Preplenish = M_Pq .* PAN_ydt; + Pejected = Pavailable.* vesicleReleaseRate(:,t)*dt; + Preprocessed = M_Pq.*Preprocess.* PAN_xdt; + + ANprobability(:,t)= min(Pejected,1); + reuptakeandlost= PAN_rdt_plus_ldt .* Pcleft; + reuptake= PAN_rdt.* Pcleft; + + Pavailable= Pavailable+ Preplenish- Pejected+ Preprocessed; + Pcleft= Pcleft + Pejected - reuptakeandlost; + Preprocess= Preprocess + reuptake - Preprocessed; + Pavailable(Pavailable<0)=0; + savePavailableSeg(:,t)=Pavailable; % synapse tracking + end + % and save it as *rate* + ANrate=ANprobability/dt; + ANprobRateOutput(:, segmentStartPTR:... + segmentStartPTR+segmentLength-1)= ANrate; + % monitor synapse contents (only sometimes used) + savePavailable(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=... + savePavailableSeg; + + % Estimate efferent effects. ARattenuation (0 <> 1) + % acoustic reflex + ARAttSeg=mean(ANrate(1:nBFs,:),1); %LSR channels are 1:nBF + % smooth + [ARAttSeg, ARboundaryProb] = ... + filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundaryProb); + ARAttSeg=ARAttSeg-ARrateThreshold; + ARAttSeg(ARAttSeg<0)=0; % prevent negative strengths + ARattenuation(segmentStartPTR:segmentEndPTR)=... + (1-ARrateToAttenuationFactorProb.* ARAttSeg); + + % MOC attenuation + % within-channel HSR response only + HSRbegins=nBFs*(nANfiberTypes-1)+1; + rates=ANrate(HSRbegins:end,:); + for idx=1:nBFs + [smoothedRates, MOCprobBoundary{idx}] = ... + filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ... + MOCprobBoundary{idx}); + smoothedRates=smoothedRates-MOCrateThreshold; + smoothedRates(smoothedRates<0)=0; + MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ... + (1- smoothedRates* rateToAttenuationFactorProb); + end + MOCattenuation(MOCattenuation<0)=0.001; + + + case 'spikes' + ANtimeCount=0; + % implement speed upt + for t = ANspeedUpFactor:ANspeedUpFactor:segmentLength; + ANtimeCount=ANtimeCount+1; + % convert release rate to probabilities + releaseProb=vesicleReleaseRate(:,t)*ANdt; + % releaseProb is the release probability per channel + % but each channel has many synapses + releaseProb=repmat(releaseProb',nFibersPerChannel,1); + releaseProb=reshape(releaseProb, nFibersPerChannel*nChannels,1); + + % AN_available=round(AN_available); % vesicles must be integer, (?needed) + M_q=AN_M- AN_available; % number of missing vesicles + M_q(M_q<0)= 0; % cannot be less than 0 + + % AN_N1 converts probability to discrete events + % it considers each event that might occur + % (how many vesicles might be released) + % and returns a count of how many were released + + % slow line +% probabilities= 1-(1-releaseProb).^AN_available; + probabilities= 1-intpow((1-releaseProb), AN_available); + ejected= probabilities> rand(length(AN_available),1); + + reuptakeandlost = AN_rdt_plus_ldt .* AN_cleft; + reuptake = AN_rdt.* AN_cleft; + + % slow line +% probabilities= 1-(1-AN_reprocess.*AN_xdt).^M_q; + probabilities= 1-intpow((1-AN_reprocess.*AN_xdt), M_q); + reprocessed= probabilities>rand(length(M_q),1); + + % slow line +% probabilities= 1-(1-AN_ydt).^M_q; + probabilities= 1-intpow((1-AN_ydt), M_q); + + replenish= probabilities>rand(length(M_q),1); + + AN_available = AN_available + replenish - ejected ... + + reprocessed; + AN_cleft = AN_cleft + ejected - reuptakeandlost; + AN_reprocess = AN_reprocess + reuptake - reprocessed; + + % ANspikes is logical record of vesicle release events>0 + ANspikes(:, ANtimeCount)= ejected; + end % t + + % zero any events that are preceded by release events ... + % within the refractory period + % The refractory period consist of two periods + % 1) the absolute period where no spikes occur + % 2) a relative period where a spike may occur. This relative + % period is realised as a variable length interval + % where the length is chosen at random + % (esentially a linear ramp up) + + % Andreas has a fix for this + for t = 1:ANtimeCount-2*lengthAbsRefractory; + % identify all spikes across fiber array at time (t) + % idx is a list of channels where spikes occurred + % ?? try sparse matrices? + idx=find(ANspikes(:,t)); + for j=idx % consider each spike + % specify variable refractory period + % between abs and 2*abs refractory period + nPointsRefractory=lengthAbsRefractory+... + round(rand*lengthAbsRefractory); + % disable spike potential for refractory period + % set all values in this range to 0 + ANspikes(j,t+1:t+nPointsRefractory)=0; + end + end %t + + % segment debugging + % plotInstructions.figureNo=98; + % plotInstructions.displaydt=ANdt; + % plotInstructions.numPlots=1; + % plotInstructions.subPlotNo=1; + % UTIL_plotMatrix(ANspikes, plotInstructions); + + % and save it. NB, AN is now on 'speedUp' time + ANoutput(:, reducedSegmentPTR: shorterSegmentEndPTR)=ANspikes; + + + %% CN Macgregor first neucleus ------------------------------- + % input is from AN so ANdt is used throughout + % Each CNneuron has a unique set of input fibers selected + % at random from the available AN fibers (CNinputfiberLists) + + % Create the dendritic current for that neuron + % First get input spikes to this neuron + synapseNo=1; + for ch=1:nChannels + for idx=1:nCNneuronsPerChannel + % determine candidate fibers for this unit + fibersUsed=CNinputfiberLists(synapseNo,:); + % ANpsth has a bin width of ANdt + % (just a simple sum across fibers) + AN_PSTH(synapseNo,:) = ... + sum(ANspikes(fibersUsed,:), 1); + synapseNo=synapseNo+1; + end + end + + % One alpha function per spike + [alphaRows alphaCols]=size(CNtrailingAlphas); + + for unitNo=1:nCNneurons + CNcurrentTemp(unitNo,:)= ... + conv(AN_PSTH(unitNo,:),CNalphaFunction); + end +% disp(['sum(AN_PSTH)= ' num2str(sum(AN_PSTH(1,:)))]) + % add post-synaptic current left over from previous segment + CNcurrentTemp(:,1:alphaCols)=... + CNcurrentTemp(:,1:alphaCols)+ CNtrailingAlphas; + + % take post-synaptic current for this segment + CNcurrentInput= CNcurrentTemp(:, 1:reducedSegmentLength); +% disp(['mean(CNcurrentInput)= ' num2str(mean(CNcurrentInput(1,:)))]) + + % trailingalphas are the ends of the alpha functions that + % spill over into the next segment + CNtrailingAlphas= ... + CNcurrentTemp(:, reducedSegmentLength+1:end); + + if CN_c>0 + % variable threshold condition (slow) + for t=1:reducedSegmentLength + CNtimeSinceLastSpike=CNtimeSinceLastSpike-ANdt; + s=CN_E>CN_Th & CNtimeSinceLastSpike<0 ; + CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike + dE =(-CN_E/CN_tauM + ... + CNcurrentInput(:,t)/CN_cap+(... + CN_Gk/CN_cap).*(CN_Ek-CN_E))*ANdt; + dGk=-CN_Gk*ANdt./tauGk + CN_b*s; + dTh=-(CN_Th-CN_Th0)*ANdt/CN_tauTh + CN_c*s; + CN_E=CN_E+dE; + CN_Gk=CN_Gk+dGk; + CN_Th=CN_Th+dTh; + CNmembranePotential(:,t)=CN_E+s.*(CN_Eb-CN_E)+CN_Er; + end + else + % static threshold (faster) + E=zeros(1,reducedSegmentLength); + Gk=zeros(1,reducedSegmentLength); + ss=zeros(1,reducedSegmentLength); + for t=1:reducedSegmentLength + % time of previous spike moves back in time + CNtimeSinceLastSpike=CNtimeSinceLastSpike-ANdt; + % action potential if E>threshold + % allow time for s to reset between events + s=CN_E>CN_Th0 & CNtimeSinceLastSpike<0 ; + ss(t)=s(1); + CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike + dE = (-CN_E/CN_tauM + ... + CNcurrentInput(:,t)/CN_cap +... + (CN_Gk/CN_cap).*(CN_Ek-CN_E))*ANdt; + dGk=-CN_Gk*ANdt./tauGk +CN_b*s; + CN_E=CN_E+dE; + CN_Gk=CN_Gk+dGk; + E(t)=CN_E(1); + Gk(t)=CN_Gk(1); + % add spike to CN_E and add resting potential (-60 mV) + CNmembranePotential(:,t)=CN_E +s.*(CN_Eb-CN_E)+CN_Er; + end + end +% disp(['CN_E= ' num2str(sum(CN_E(1,:)))]) +% disp(['CN_Gk= ' num2str(sum(CN_Gk(1,:)))]) +% disp(['CNmembranePotential= ' num2str(sum(CNmembranePotential(1,:)))]) +% plot(CNmembranePotential(1,:)) + + + % extract spikes. A spike is a substantial upswing in voltage + CN_spikes=CNmembranePotential> -0.02; +% disp(['CNspikesbefore= ' num2str(sum(sum(CN_spikes)))]) + + % now remove any spike that is immediately followed by a spike + % NB 'find' works on columns (whence the transposing) + % for each spike put a zero in the next epoch + CN_spikes=CN_spikes'; + idx=find(CN_spikes); + idx=idx(1:end-1); + CN_spikes(idx+1)=0; + CN_spikes=CN_spikes'; +% disp(['CNspikes= ' num2str(sum(sum(CN_spikes)))]) + + % segment debugging + % plotInstructions.figureNo=98; + % plotInstructions.displaydt=ANdt; + % plotInstructions.numPlots=1; + % plotInstructions.subPlotNo=1; + % UTIL_plotMatrix(CN_spikes, plotInstructions); + + % and save it + CNoutput(:, reducedSegmentPTR:shorterSegmentEndPTR)=... + CN_spikes; + + + %% IC ---------------------------------------------- + % MacGregor or some other second order neurons + + % combine CN neurons in same channel, i.e. same BF & same tauCa + % to generate inputs to single IC unit + channelNo=0; + for idx=1:nCNneuronsPerChannel:nCNneurons-nCNneuronsPerChannel+1; + channelNo=channelNo+1; + CN_PSTH(channelNo,:)=... + sum(CN_spikes(idx:idx+nCNneuronsPerChannel-1,:)); + end + + [alphaRows alphaCols]=size(ICtrailingAlphas); + for ICneuronNo=1:nICcells + ICcurrentTemp(ICneuronNo,:)= ... + conv(CN_PSTH(ICneuronNo,:), IC_CNalphaFunction); + end + + % add the unused current from the previous convolution + ICcurrentTemp(:,1:alphaCols)=ICcurrentTemp(:,1:alphaCols)... + + ICtrailingAlphas; + % take what is required and keep the trailing part for next time + inputCurrent=ICcurrentTemp(:, 1:reducedSegmentLength); + ICtrailingAlphas=ICcurrentTemp(:, reducedSegmentLength+1:end); + + if IC_c==0 + % faster computation when threshold is stable (C==0) + for t=1:reducedSegmentLength + s=IC_E>IC_Th0; + dE = (-IC_E/IC_tauM + inputCurrent(:,t)/IC_cap +... + (IC_Gk/IC_cap).*(IC_Ek-IC_E))*ANdt; + dGk=-IC_Gk*ANdt/IC_tauGk +IC_b*s; + IC_E=IC_E+dE; + IC_Gk=IC_Gk+dGk; + ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er; + end + else + % threshold is changing (IC_c>0; e.g. bushy cell) + for t=1:reducedSegmentLength + dE = (-IC_E/IC_tauM + ... + inputCurrent(:,t)/IC_cap + (IC_Gk/IC_cap)... + .*(IC_Ek-IC_E))*ANdt; + IC_E=IC_E+dE; + s=IC_E>IC_Th; + ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er; + dGk=-IC_Gk*ANdt/IC_tauGk +IC_b*s; + IC_Gk=IC_Gk+dGk; + + % After a spike, the threshold is raised + % otherwise it settles to its baseline + dTh=-(IC_Th-Th0)*ANdt/IC_tauTh +s*IC_c; + IC_Th=IC_Th+dTh; + end + end + + ICspikes=ICmembranePotential> -0.01; + % now remove any spike that is immediately followed by a spike + % NB 'find' works on columns (whence the transposing) + ICspikes=ICspikes'; + idx=find(ICspikes); + idx=idx(1:end-1); + ICspikes(idx+1)=0; + ICspikes=ICspikes'; + + nCellsPerTau= nICcells/nANfiberTypes; + firstCell=1; + lastCell=nCellsPerTau; + for tauCount=1:nANfiberTypes + % separate rates according to fiber types + % currently only the last segment is saved + ICfiberTypeRates(tauCount, ... + reducedSegmentPTR:shorterSegmentEndPTR)=... + sum(ICspikes(firstCell:lastCell, :))... + /(nCellsPerTau*ANdt); + firstCell=firstCell+nCellsPerTau; + lastCell=lastCell+nCellsPerTau; + end + + ICoutput(:,reducedSegmentPTR:shorterSegmentEndPTR)=ICspikes; + + % store membrane output on original dt scale + if nBFs==1 % single channel + x= repmat(ICmembranePotential(1,:), ANspeedUpFactor,1); + x= reshape(x,1,segmentLength); + if nANfiberTypes>1 % save HSR and LSR + y=repmat(ICmembranePotential(end,:),... + ANspeedUpFactor,1); + y= reshape(y,1,segmentLength); + x=[x; y]; + end + ICmembraneOutput(:, segmentStartPTR:segmentEndPTR)= x; + end + + % estimate efferent effects. + % ARis based on LSR units. LSR channels are 1:nBF + if nANfiberTypes>1 % AR is multi-channel only + ARAttSeg=sum(ICspikes(1:nBFs,:),1)/ANdt; + [ARAttSeg, ARboundary] = ... + filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundary); + ARAttSeg=ARAttSeg-ARrateThreshold; + ARAttSeg(ARAttSeg<0)=0; % prevent negative strengths + % scale up to dt from ANdt + x= repmat(ARAttSeg, ANspeedUpFactor,1); + x=reshape(x,1,segmentLength); + ARattenuation(segmentStartPTR:segmentEndPTR)=... + (1-ARrateToAttenuationFactor* x); + ARattenuation(ARattenuation<0)=0.001; + else + % single channel model; disable AR + ARattenuation(segmentStartPTR:segmentEndPTR)=... + ones(1,segmentLength); + end + + % MOC attenuation using HSR response only + % Separate MOC effect for each BF + HSRbegins=nBFs*(nANfiberTypes-1)+1; + rates=ICspikes(HSRbegins:end,:)/ANdt; + for idx=1:nBFs + [smoothedRates, MOCboundary{idx}] = ... + filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ... + MOCboundary{idx}); + MOCattSegment(idx,:)=smoothedRates; + % expand timescale back to model dt from ANdt + x= repmat(MOCattSegment(idx,:), ANspeedUpFactor,1); + x= reshape(x,1,segmentLength); + MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ... + (1- MOCrateToAttenuationFactor* x); + end + MOCattenuation(MOCattenuation<0)=0.04; + % segment debugging + % plotInstructions.figureNo=98; + % plotInstructions.displaydt=ANdt; + % plotInstructions.numPlots=1; + % plotInstructions.subPlotNo=1; + % UTIL_plotMatrix(ICspikes, plotInstructions); + + end % AN_spikesOrProbability + segmentStartPTR=segmentStartPTR+segmentLength; + reducedSegmentPTR=reducedSegmentPTR+reducedSegmentLength; + + +end % segment + +path(restorePath)
--- a/parameterStore/MAPparamsEndo.m Mon Jun 06 09:11:29 2011 +0100 +++ b/parameterStore/MAPparamsEndo.m Tue Jun 07 09:53:50 2011 +0100 @@ -137,7 +137,7 @@ IHC_cilia_RPParams.Cab= 4e-012; % IHC capacitance (F) IHC_cilia_RPParams.Cab= 1e-012; % IHC capacitance (F) IHC_cilia_RPParams.Et= 0.100; % endocochlear potential (V) -IHC_cilia_RPParams.Et= 0.07; % endocochlear potential (V) +IHC_cilia_RPParams.Et= 0.065; % 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
--- a/parameterStore/MAPparamsNormal.m Mon Jun 06 09:11:29 2011 +0100 +++ b/parameterStore/MAPparamsNormal.m Tue Jun 07 09:53:50 2011 +0100 @@ -20,7 +20,6 @@ currentFile=mfilename; % i.e. the name of this mfile method.parameterSource=currentFile(10:end); % for the record -switchOffEfferent=0; efferentDelay=0.010; method.segmentDuration=efferentDelay; @@ -57,17 +56,14 @@ % Acoustic reflex: maximum attenuation should be around 25 dB Price (1966) % i.e. a minimum ratio of 0.056. -if ~switchOffEfferent - % 'spikes' model: AR based on brainstem spiking activity (LSR) - OMEParams.rateToAttenuationFactor=0.003; % * N(all ICspikes) +% 'spikes' model: AR based on brainstem spiking activity (LSR) +OMEParams.rateToAttenuationFactor=0.003; % * N(all ICspikes) % OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) - % 'probability model': Ar based on An firing probabilities (LSR) - OMEParams.rateToAttenuationFactorProb=0.005;% * N(all ANrates) + +% 'probability model': Ar based on AN firing probabilities (LSR) +OMEParams.rateToAttenuationFactorProb=0.003;% * N(all ANrates) % OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates) -else - OMEParams.rateToAttenuationFactor=0; % 0= off - OMEParams.rateToAttenuationFactorProb=0; % 0= off -end + % asymptote should be around 100-200 ms OMEParams.ARtau=.05; % AR smoothing function % delay must be longer than the segment length @@ -104,18 +100,14 @@ % DRNL MOC efferents DRNLParams.MOCdelay = efferentDelay; % must be < segment length! -if ~switchOffEfferent - % 'spikes' model: MOC based on brainstem spiking activity (HSR) - DRNLParams.rateToAttenuationFactor = .009; % strength of MOC - DRNLParams.rateToAttenuationFactor = .009; % strength of MOC +% 'spikes' model: MOC based on brainstem spiking activity (HSR) +DRNLParams.rateToAttenuationFactor = .009; % strength of MOC +DRNLParams.rateToAttenuationFactor = .009; % strength of MOC % DRNLParams.rateToAttenuationFactor = 0; % strength of MOC - % 'spikes' model: MOC based on brainstem spiking activity (HSR) - DRNLParams.rateToAttenuationFactorProb = .002; % strength of MOC -else - DRNLParams.rateToAttenuationFactor = 0; % 0 = MOC off (probability) - DRNLParams.rateToAttenuationFactorProb = 0; % 0 = MOC off (spikes) -end +% 'probability' model: MOC based on AN spiking activity (HSR) +DRNLParams.rateToAttenuationFactorProb = .007; % strength of MOC +DRNLParams.rateToAttenuationFactorProb = .0; % strength of MOC DRNLParams.MOCtau =.03; % smoothing for MOC DRNLParams.MOCrateThreshold =50; % set to AN rate threshold @@ -137,7 +129,6 @@ IHC_cilia_RPParams.Cab= 4e-012; % IHC capacitance (F) IHC_cilia_RPParams.Cab= 1e-012; % IHC capacitance (F) IHC_cilia_RPParams.Et= 0.100; % endocochlear potential (V) -IHC_cilia_RPParams.Et= 0.07; % 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
--- a/parameterStore/MAPparamsOHC.m Mon Jun 06 09:11:29 2011 +0100 +++ b/parameterStore/MAPparamsOHC.m Tue Jun 07 09:53:50 2011 +0100 @@ -20,7 +20,6 @@ currentFile=mfilename; % i.e. the name of this mfile method.parameterSource=currentFile(10:end); % for the record -switchOffEfferent=0; efferentDelay=0.010; method.segmentDuration=efferentDelay; @@ -57,17 +56,14 @@ % Acoustic reflex: maximum attenuation should be around 25 dB Price (1966) % i.e. a minimum ratio of 0.056. -if ~switchOffEfferent - % 'spikes' model: AR based on brainstem spiking activity (LSR) - OMEParams.rateToAttenuationFactor=0.003; % * N(all ICspikes) +% 'spikes' model: AR based on brainstem spiking activity (LSR) +OMEParams.rateToAttenuationFactor=0.003; % * N(all ICspikes) % OMEParams.rateToAttenuationFactor=0; % * N(all ICspikes) - % 'probability model': Ar based on An firing probabilities (LSR) - OMEParams.rateToAttenuationFactorProb=0.005;% * N(all ANrates) + +% 'probability model': Ar based on AN firing probabilities (LSR) +OMEParams.rateToAttenuationFactorProb=0.005;% * N(all ANrates) % OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates) -else - OMEParams.rateToAttenuationFactor=0; % 0= off - OMEParams.rateToAttenuationFactorProb=0; % 0= off -end + % asymptote should be around 100-200 ms OMEParams.ARtau=.05; % AR smoothing function % delay must be longer than the segment length @@ -104,18 +100,13 @@ % DRNL MOC efferents DRNLParams.MOCdelay = efferentDelay; % must be < segment length! -if ~switchOffEfferent - % 'spikes' model: MOC based on brainstem spiking activity (HSR) - DRNLParams.rateToAttenuationFactor = .009; % strength of MOC - DRNLParams.rateToAttenuationFactor = .009; % strength of MOC +% 'spikes' model: MOC based on brainstem spiking activity (HSR) +DRNLParams.rateToAttenuationFactor = .009; % strength of MOC +DRNLParams.rateToAttenuationFactor = .009; % strength of MOC % DRNLParams.rateToAttenuationFactor = 0; % strength of MOC - % 'spikes' model: MOC based on brainstem spiking activity (HSR) - DRNLParams.rateToAttenuationFactorProb = .002; % strength of MOC -else - DRNLParams.rateToAttenuationFactor = 0; % 0 = MOC off (probability) - DRNLParams.rateToAttenuationFactorProb = 0; % 0 = MOC off (spikes) -end +% 'probability' model: MOC based on AN spiking activity (HSR) +DRNLParams.rateToAttenuationFactorProb = .002; % strength of MOC DRNLParams.MOCtau =.03; % smoothing for MOC DRNLParams.MOCrateThreshold =50; % set to AN rate threshold @@ -137,7 +128,6 @@ IHC_cilia_RPParams.Cab= 4e-012; % IHC capacitance (F) IHC_cilia_RPParams.Cab= 1e-012; % IHC capacitance (F) IHC_cilia_RPParams.Et= 0.100; % endocochlear potential (V) -% IHC_cilia_RPParams.Et= 0.07; % 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 @@ -183,6 +173,8 @@ 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 @@ -212,13 +204,13 @@ MacGregorMultiParams.dendriteLPfreq=50; % dendritic filter MacGregorMultiParams.currentPerSpike=35e-9; % *per spike -% MacGregorMultiParams.currentPerSpike=45e-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.0001;% K+ conductance tau (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 @@ -233,12 +225,13 @@ 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.0003; % K+ conductance tau (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
--- a/testPrograms/showMAP.m Mon Jun 06 09:11:29 2011 +0100 +++ b/testPrograms/showMAP.m Tue Jun 07 09:53:50 2011 +0100 @@ -1,10 +1,12 @@ function showMAP (options) % defaults -% options.showModelParameters=1; -% options.showModelOutput=1; -% options.printFiringRates=1; -% options.showACF=1; -% options.showEfferent=1; +% options.showModelParameters=1; +% options.showModelOutput=1; +% options.printFiringRates=1; +% options.showACF=1; +% options.showEfferent=1; +% options.surfProbability=0; +% options.fileName=[]; dbstop if warning @@ -13,6 +15,9 @@ DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV... IHCoutput ANprobRateOutput ANoutput savePavailable tauCas ... CNoutput ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation +global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams +global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams + restorePath=path; addpath ( ['..' filesep 'utilities'], ['..' filesep 'parameterStore']) @@ -21,8 +26,10 @@ options.showModelParameters=1; options.showModelOutput=1; options.printFiringRates=1; - options.showACF=1; + options.showACF=0; options.showEfferent=1; + options.surfProbability=0; + options.fileName=[]; end if options.showModelParameters @@ -46,13 +53,13 @@ duration=size(TMoutput,2)*dt; disp(['AN: ' num2str(sum(sum(ANoutput(end-nHSRfibers+1:end,:)))/... (nHSRfibers*duration))]) - + nCNneurons=size(CNoutput,1); nHSRCNneuronss=nCNneurons/nANfiberTypes; disp(['CN: ' num2str(sum(sum(CNoutput(end-nHSRCNneuronss+1:end,:)))... /(nHSRCNneuronss*duration))]) disp(['IC: ' num2str(sum(sum(ICoutput))/duration)]) -% disp(['IC by type: ' num2str(mean(ICfiberTypeRates,2)')]) + % disp(['IC by type: ' num2str(mean(ICfiberTypeRates,2)')]) else disp(['AN: ' num2str(mean(mean(ANprobRateOutput)))]) end @@ -64,7 +71,7 @@ plotInstructions.figureNo=99; signalRMS=mean(savedInputSignal.^2)^0.5; signalRMSdb=20*log10(signalRMS/20e-6); - + % plot signal (1) plotInstructions.displaydt=dt; plotInstructions.numPlots=6; @@ -74,18 +81,18 @@ r=size(savedInputSignal,1); if r==1, savedInputSignal=savedInputSignal'; end UTIL_plotMatrix(savedInputSignal', plotInstructions); - + % stapes (2) plotInstructions.subPlotNo=2; plotInstructions.title= ['stapes displacement']; UTIL_plotMatrix(OMEoutput, plotInstructions); - + % DRNL (3) plotInstructions.subPlotNo=3; plotInstructions.title= ['BM displacement']; plotInstructions.yValues= savedBFlist; UTIL_plotMatrix(DRNLoutput, plotInstructions); - + switch saveAN_spikesOrProbability case 'spikes' % AN (4) @@ -100,7 +107,7 @@ plotInstructions.rasterDotSize=3; end UTIL_plotMatrix(ANoutput, plotInstructions); - + % CN (5) plotInstructions.displaydt=ANdt; plotInstructions.subPlotNo=5; @@ -109,7 +116,7 @@ plotInstructions.rasterDotSize=3; end UTIL_plotMatrix(CNoutput, plotInstructions); - + % IC (6) plotInstructions.displaydt=ANdt; plotInstructions.subPlotNo=6; @@ -126,7 +133,7 @@ plotInstructions.zValuesRange= [-.1 0]; UTIL_plotMatrix(ICmembraneOutput, plotInstructions); end - + otherwise % probability (4-6) plotInstructions.displaydt=dt; plotInstructions.numPlots=2; @@ -153,7 +160,7 @@ plotInstructions.title= ['AR strength. Signal level= ' ... num2str(signalRMSdb,'%4.0f') ' dB SPL']; UTIL_plotMatrix(20*log10(ARattenuation), plotInstructions); - + plotInstructions.subPlotNo=2; plotInstructions.yValues= savedBFlist; plotInstructions.yLabel= 'BF'; @@ -165,6 +172,24 @@ colorbar end +%% surface plot of probability +if options.surfProbability + figure(97), clf + % select only HSR fibers at the bottom of the matrix + ANprobRateOutput= ANprobRateOutput(end-length(savedBFlist)+1:end,:); + [nY nX]=size(ANprobRateOutput); + time=dt*(1:nX); + surf(time, savedBFlist, ANprobRateOutput) + shading interp + set(gca, 'yScale','log') + xlim([0 max(time)]), ylim([0 max(savedBFlist)]), zlim([0 1000]) + xlabel('time (s)') + ylabel('best frequency (Hz)') + zlabel('spike rate') + view([-20 60]) + title ([options.fileName ': ' num2str(signalRMSdb,'% 3.0f') ' dB']) +end + %% ACF plot if required if options.showACF @@ -172,7 +197,7 @@ method.dt=dt; method.segmentNo=1; method.nonlinCF=savedBFlist; - + minPitch= 80; maxPitch= 4000; numPitches=100; % specify lags pitches=10.^ linspace(log10(minPitch), log10(maxPitch),numPitches); pitches=fliplr(pitches); @@ -180,18 +205,18 @@ filteredSACFParams.acfTau= .003; % time constant of running ACF filteredSACFParams.lambda= 0.12; % slower filter to smooth ACF filteredSACFParams.lambda= 0.01; % slower filter to smooth ACF - + filteredSACFParams.plotACFs=0; % special plot (see code) filteredSACFParams.plotFilteredSACF=0; % 0 plots unfiltered ACFs filteredSACFParams.plotMoviePauses=.3; % special plot (see code) - + filteredSACFParams.usePressnitzer=0; % attenuates ACF at long lags filteredSACFParams.lagsProcedure= 'useAllLags'; % filteredSACFParams.lagsProcedure= 'useBernsteinLagWeights'; % filteredSACFParams.lagsProcedure= 'omitShortLags'; filteredSACFParams.criterionForOmittingLags=3; filteredSACFParams.plotACFsInterval=200; - + if filteredSACFParams.plotACFs % plot original waveform on ACF plot figure(13), clf @@ -201,16 +226,16 @@ xlim([0 t(end)]) title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); end - + % plot original waveform on summary/smoothed ACF plot - figure(97), clf + figure(96), clf subplot(2,1,1) t=dt*(1:length(savedInputSignal)); plot(t,savedInputSignal) xlim([0 t(end)]) title(['stimulus: ' num2str(signalRMSdb, '%4.0f') ' dB SPL']); - - + + % compute ACF switch saveAN_spikesOrProbability case 'probability' @@ -218,12 +243,12 @@ otherwise inputToACF=ANoutput; end - + disp ('computing ACF...') [P, BFlist, sacf, boundaryValue] = ... filteredSACF(inputToACF, method, filteredSACFParams); disp(' ACF done.') - + % SACF subplot(2,1,2) imagesc(P)
--- a/testPrograms/test_MAP1_14.m Mon Jun 06 09:11:29 2011 +0100 +++ b/testPrograms/test_MAP1_14.m Tue Jun 07 09:53:50 2011 +0100 @@ -40,7 +40,7 @@ %% #2 probability (fast) or spikes (slow) representation AN_spikesOrProbability='spikes'; % or -% AN_spikesOrProbability='probability'; +AN_spikesOrProbability='probability'; %% #3 pure tone, harmonic sequence or speech file input @@ -54,14 +54,19 @@ rampDuration=.005; % seconds % or -% signalType= 'file'; -% fileName='twister_44kHz'; +signalType= 'file'; +fileName='twister_44kHz'; % fileName='new-da-44khz'; +% mix with an optional second file? +mixerFile=[]; +%or +mixerFile='babble'; +leveldBSPL2=-60; %% #4 rms level % signal details -leveldBSPL= 90; % dB SPL +leveldBSPL= 60; % dB SPL %% #5 number of channels in the model @@ -71,7 +76,7 @@ BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels)); % or specify your own channel BFs -BFlist=toneFrequency; +% BFlist=toneFrequency; %% #6 change model parameters @@ -82,18 +87,26 @@ % *after* the MAPparams file has been read % This example declares only one fiber type with a calcium clearance time % constant of 80e-6 s (HSR fiber) when the probability option is selected. - paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ... - 'IHCpreSynapseParams.tauCa=86e-6;'}; +% paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ... +% 'IHCpreSynapseParams.tauCa=86e-6;'}; +% paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ... +% 'DRNLParams.rateToAttenuationFactorProb = 0;'}; %% delare showMap options showMapOptions=[]; % use defaults % or (example: show everything including an smoothed SACF output - showMapOptions.showModelParameters=0; - showMapOptions.showModelOutput=1; - showMapOptions.printFiringRates=1; - showMapOptions.showACF=0; - showMapOptions.showEfferent=1; +showMapOptions.showModelParameters=1; +showMapOptions.showModelOutput=1; +showMapOptions.printFiringRates=1; +showMapOptions.showACF=0; +showMapOptions.showEfferent=1; +if strcmp(AN_spikesOrProbability, 'probability') + showMapOptions.surfProbability=1; +end +if strcmp(signalType, 'file') + showMapOptions.fileName=fileName; +end %% Generate stimuli @@ -106,12 +119,24 @@ leveldBSPL, duration, rampDuration); case 'file' + %% file input simple or mixed [inputSignal sampleRate]=wavread(fileName); - inputSignal(:,1); + inputSignal=inputSignal(:,1); targetRMS=20e-6*10^(leveldBSPL/20); rms=(mean(inputSignal.^2))^0.5; amp=targetRMS/rms; inputSignal=inputSignal*amp; + if ~isempty(mixerFile) + [inputSignal2 sampleRate]=wavread(mixerFile); + inputSignal2=inputSignal2(:,1); + [r c]=size(inputSignal); + inputSignal2=inputSignal2(1:r); + targetRMS=20e-6*10^(leveldBSPL2/20); + rms=(mean(inputSignal2.^2))^0.5; + amp=targetRMS/rms; + inputSignal2=inputSignal2*amp; + inputSignal=inputSignal+inputSignal2; + end end