Mercurial > hg > map
diff userProgramsTim/MAP1_14_olde_version_dont_use.m @ 38:c2204b18f4a2 tip
End nov big change
| author | Ray Meddis <rmeddis@essex.ac.uk> |
|---|---|
| date | Mon, 28 Nov 2011 13:34:28 +0000 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/userProgramsTim/MAP1_14_olde_version_dont_use.m Mon Nov 28 13:34:28 2011 +0000 @@ -0,0 +1,1269 @@ +function MAP1_14(inputSignal, sampleRate, BFlist, MAPparamsName, ... + AN_spikesOrProbability, paramChanges) +% To test this function use test_MAP1_14 in this folder +% +% example: +% <navigate to 'MAP1_14\MAP'> +% [inputSignal FS] = wavread('../wavFileStore/twister_44kHz'); +% MAP1_14(inputSignal, FS, -1, 'Normal', 'probability', []) +% +% All arguments are mandatory. +% +% BFlist is a vector 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 +% e.g. paramChanges{1}= 'DRNLParams.a=0;'; % disable OHCs +% e.g. paramchanges={}; % no changes +% 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 ANtauCas ... + CNtauGk 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. + +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 +% Beware, 'BFlist=-1' is a legitimate argument for MAPparams<> +% It means that the calling program allows MAPparams to specify the list +cmd=['method=MAPparams' MAPparamsName ... + '(BFlist, sampleRate, 0, paramChanges);']; +eval(cmd); +BFlist=DRNLParams.nonlinCFs; + +% save as global for later plotting if required +savedBFlist=BFlist; +saveAN_spikesOrProbability=AN_spikesOrProbability; +saveMAPparamsName=MAPparamsName; + +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; +MOCrateThresholdProb=DRNLParams.MOCrateThresholdProb; + +% 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 +DRNLcompressionThreshold=DRNLParams.cTh; +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); +ANtauCas= IHCpreSynapseParams.tauCa; +nANchannels= nANfiberTypes*nBFs; +synapticCa= zeros(nANchannels,segmentLength); + +% Calcium control (more calcium, greater release rate) +ECa=IHCpreSynapseParams.ECa; +gamma=IHCpreSynapseParams.gamma; +beta=IHCpreSynapseParams.beta; +tauM=IHCpreSynapseParams.tauM; +mICa=zeros(nANchannels,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 ANtauCas (one pre fiber type) +tauCa=repmat(ANtauCas, nBFs,1); +tauCa=reshape(tauCa, nANchannels, 1); + +% presynapse startup values (vectors, length:nANchannels) +% 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= nANchannels*nFibersPerChannel; +AN_refractory_period= AN_IHCsynapseParams.refractory_period; + +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, nANchannels,1); +PAN_ldt = repmat(AN_IHCsynapseParams.l*dt, nANchannels,1); +PAN_xdt = repmat(AN_IHCsynapseParams.x*dt, nANchannels,1); +PAN_rdt = repmat(AN_IHCsynapseParams.r*dt, nANchannels,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(nANchannels,segmentLength); +ANprobRateOutput=zeros(nANchannels,signalLength); +lengthAbsRefractoryP= round(AN_refractory_period/dt); +cumANnotFireProb=ones(nANchannels,signalLength); +% special variables for monitoring synaptic cleft (specialists only) +savePavailableSeg=zeros(nANchannels,segmentLength); +savePavailable=zeros(nANchannels,signalLength); + +% spikes % ! ! ! ! ! ! ! ! +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 +% CNtauGk (Potassium time constant) determines the rate of firing of +% the unit when driven hard by a DC input (not normally >350 sp/s) +% If there is more than one value, everything is replicated accordingly + +ANavailableFibersPerChan=AN_IHCsynapseParams.numFibers; +ANfibersFanInToCN=MacGregorMultiParams.fibersPerNeuron; + +CNtauGk=MacGregorMultiParams.tauGk; % row vector of CN types (by tauGk) +nCNtauGk=length(CNtauGk); + +% the total number of 'channels' is now greater +% 'channel' is defined as collections of units with the same parameters +% i.e. same BF, same ANtau, same CNtauGk +nCNchannels=nANchannels*nCNtauGk; + +nCNneuronsPerChannel=MacGregorMultiParams.nNeuronsPerBF; +tauGk=repmat(CNtauGk, nCNneuronsPerChannel,1); +tauGk=reshape(tauGk,nCNneuronsPerChannel*nCNtauGk,1); + +% Now the number of neurons has been increased +nCNneurons=nCNneuronsPerChannel*nCNchannels; +CNmembranePotential=zeros(nCNneurons,reducedSegmentLength); + +% establish which ANfibers (by name) feed into which CN nuerons +CNinputfiberLists=zeros(nANchannels*nCNneuronsPerChannel, ANfibersFanInToCN); +unitNo=1; +for ch=1:nANchannels + % Each channel contains a number of units =length(listOfFanInValues) + for idx=1:nCNneuronsPerChannel + for idx2=1:nCNtauGk + fibersUsed=(ch-1)*ANavailableFibersPerChan + ... + ceil(rand(1,ANfibersFanInToCN)* ANavailableFibersPerChan); + CNinputfiberLists(unitNo,:)=fibersUsed; + unitNo=unitNo+1; + end + 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(tauGk,nANchannels,1); + +CNoutput=false(nCNneurons,reducedSignalLength); + + +%% MacGregor (IC - second nucleus) -------- +nICcells=nANchannels*nCNtauGk; % one cell per channel +CN_PSTH=zeros(nICcells ,reducedSegmentLength); + +% 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(nICcells,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; + + inputPressureSegment=inputSignal... + (:,segmentStartPTR:segmentStartPTR+segmentLength-1); + + % segment debugging plots + % figure(98) + % plot(segmentTime,inputPressureSegment), title('signalSegment') + + + % OME ---------------------- + + % OME Stage 1: external resonances. Add to inputSignal pressure wave + y=inputPressureSegment; + for n=1:nOMEExtFilters + % any number of resonances can be used + [x OMEExtFilterBndry{n}] = ... + filter(ExtFilter_b{n},ExtFilter_a{n},... + inputPressureSegment, OMEExtFilterBndry{n}); + x= x* OMEgainScalars(n); + % This is a parallel resonance so add it + y=y+x; + end + inputPressureSegment=y; + OMEextEarPressure(segmentStartPTR:segmentEndPTR)= inputPressureSegment; + + % OME stage 2: convert input pressure (velocity) to + % tympanic membrane(TM) displacement using low pass filter + [TMdisplacementSegment OME_TMdisplacementBndry] = ... + filter(TMdisp_b,TMdisp_a,inputPressureSegment, ... + 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 + % apply MOC to nonlinear input function + nonlinOutput=stapesDisplacement.* MOC; + + % first gammatone filter (nonlin path) + for order = 1 : GTnonlinOrder + [nonlinOutput GTnonlinBdry1{BFno,order}] = ... + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ... + nonlinOutput, GTnonlinBdry1{BFno,order}); + end + +% % original broken stick instantaneous compression +% y= nonlinOutput.* DRNLa; % linear section. +% % compress parts of the signal above the compression threshold +% abs_x = abs(nonlinOutput); +% idx=find(abs_x>DRNLcompressionThreshold); +% if ~isempty(idx)>0 +% y(idx)=sign(y(idx)).* (DRNLb*abs_x(idx).^DRNLc); +% end +% nonlinOutput=y; + + + % new broken stick instantaneous compression + y= nonlinOutput.* DRNLa; % linear section attenuation/gain. + % compress parts of the signal above the compression threshold +% holdY=y; + abs_y = abs(y); + idx=find(abs_y>DRNLcompressionThreshold); + if ~isempty(idx)>0 +% y(idx)=sign(y(idx)).* (DRNLcompressionThreshold + ... +% (abs_y(idx)-DRNLcompressionThreshold).^DRNLc); + y(idx)=sign(y(idx)).* (DRNLcompressionThreshold + ... + (abs_y(idx)-DRNLcompressionThreshold)*DRNLc); + end + nonlinOutput=y; + +% % Boltzmann compression +% y=(nonlinOutput * DRNLa*100000); +% holdY=y; +% y=abs(y); +% s=10; u=0.0; +% x=1./(1+exp(-(y-u)/s))-0.5; +% nonlinOutput=sign(nonlinOutput).*x/10000; + + +% if segmentStartPTR==10*segmentLength+1 +% figure(90) +% plot(holdY,'b'), hold on +% plot(nonlinOutput, 'r'), hold off +% ylim([-1e-5 1e-5]) +% pause(1) +% end + + % 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; +% disp(num2str(max(linOutput))) + 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 synapticCa 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; + + %% Apply refractory effect + % the probability of a spike's occurring in the preceding + % refractory window (t= tnow-refractory period to tnow-) + % pFired= 1 - II(1-p(t)), + % we need a running account of cumProb=II(1-p(t)) + % cumProb(t)= cumProb(t-1)*(1-p(t))/(1-p(t-refracPeriod)) + % cumProb(0)=0 + % pFired(t)= 1-cumProb(t) + % This gives the fraction of firing events that must be + % discounted because of a firing event in the refractory + % period + % p(t)= ANprobOutput(t) * pFired(t) + % where ANprobOutput is the uncorrected firing probability + % based on vesicle release rate + % NB this covers only the absoute refractory period + % not the relative refractory period. To approximate this it + % is necessary to extend the refractory period by 50% + + + for t = segmentStartPTR:segmentEndPTR; + if t>1 + ANprobRateOutput(:,t)= ANprobRateOutput(:,t)... + .* cumANnotFireProb(:,t-1); + end + % add recent and remove distant probabilities + refrac=round(lengthAbsRefractoryP * 1.5); + if t>refrac + cumANnotFireProb(:,t)= cumANnotFireProb(:,t-1)... + .*(1-ANprobRateOutput(:,t)*dt)... + ./(1-ANprobRateOutput(:,t-refrac)*dt); + end + end +% figure(88), plot(cumANnotFireProb'), title('cumNotFire') +% figure(89), plot(ANprobRateOutput'), title('ANprobRateOutput') + + %% Estimate efferent effects. ARattenuation (0 <> 1) + % acoustic reflex + [r c]=size(ANrate); + if r>nBFs % Only if LSR fibers are computed + 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); + end + % plot(ARattenuation) + + % MOC attenuation based on within-channel HSR fiber activity + HSRbegins=nBFs*(nANfiberTypes-1)+1; + rates=ANrate(HSRbegins:end,:); + if rateToAttenuationFactorProb<0 + % negative factor implies a fixed attenuation + MOCattenuation(:,segmentStartPTR:segmentEndPTR)= ... + ones(size(rates))* -rateToAttenuationFactorProb; + else + for idx=1:nBFs + [smoothedRates, MOCprobBoundary{idx}] = ... + filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ... + MOCprobBoundary{idx}); + smoothedRates=smoothedRates-MOCrateThresholdProb; + smoothedRates=max(smoothedRates, 0); + + x=(1- smoothedRates* rateToAttenuationFactorProb); + x=max(x, 10^(-30/20)); + MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= x; + end + end + MOCattenuation(MOCattenuation<0)=0.001; + + % plot(MOCattenuation) + + + 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*nANchannels,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:nCNchannels + 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,:)= ... + conv2(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 (BF x AN tau x CNtau) + % i.e. same BF, same tauCa, same CNtau + % 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,:)= ... + conv2(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 + % do this for single channel models only + % and only for the HSR-driven IC cells + if round(nICcells/length(ANtauCas))==1 % single channel + % select HSR driven cells + x= ICmembranePotential(length(ANtauCas),:); + % restore original dt + x= repmat(x, 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}); + % spont 'rates' is zero for IC + 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 + +%% apply refractory correction to spike probabilities + +% the probability of a spike's having occurred in the preceding +% refractory window +% pFired= 1 - II(1-p(t)), t= tnow-refractory period: tnow-1 +% we need a running account of cumProb=II(1-p(t)) +% cumProb(t)= cumProb(t-1)*(1-p(t-1))/(1-p(t-refracPeriod)) +% cumProb(0)=0 +% pFired(t)= 1-cumProb(t) +% whence +% p(t)= ANprobOutput(t) * pFired(t) +% where ANprobOutput is the uncorrected probability + + +% switch AN_spikesOrProbability +% case 'probability' +% ANprobOutput=ANprobRateOutput*dt; +% [r nEpochs]=size(ANprobOutput); +% % find probability of no spikes in refractory period +% pNoSpikesInRefrac=ones(size(ANprobOutput)); +% pSpike=zeros(size(ANprobOutput)); +% for epochNo=lengthAbsRefractoryP+2:nEpochs +% 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)