rmeddis@38: function  MAP1_14(inputSignal, sampleRate, BFlist, MAPparamsName, ...
rmeddis@38:     AN_spikesOrProbability, paramChanges)
rmeddis@38: % To test this function use test_MAP1_14 in this folder
rmeddis@38: %
rmeddis@38: % example:
rmeddis@38: % <navigate to 'MAP1_14\MAP'>
rmeddis@38: %  [inputSignal FS] = wavread('../wavFileStore/twister_44kHz');
rmeddis@38: %  MAP1_14(inputSignal, FS, -1, 'Normal', 'probability', [])
rmeddis@38: %
rmeddis@38: % All arguments are mandatory.
rmeddis@38: %
rmeddis@38: %  BFlist is a vector of BFs but can be '-1' to allow MAPparams to choose
rmeddis@38: %  MAPparamsName='Normal';          % source of model parameters
rmeddis@38: %  AN_spikesOrProbability='spikes'; % or 'probability'
rmeddis@38: %  paramChanges is a cell array of strings that can be used to make last
rmeddis@38: %   minute parameter changes, e.g., to simulate OHC loss
rmeddis@38: %  e.g.  paramChanges{1}= 'DRNLParams.a=0;'; % disable OHCs
rmeddis@38: %  e.g.  paramchanges={};                    % no changes
rmeddis@38: % The model parameters are established in the MAPparams<***> file
rmeddis@38: %  and stored as global
rmeddis@38: 
rmeddis@38: restorePath=path;
rmeddis@38: addpath (['..' filesep 'parameterStore'])
rmeddis@38: 
rmeddis@38: global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams
rmeddis@38: global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams
rmeddis@38: 
rmeddis@38: % All of the results of this function are stored as global
rmeddis@38: global dt ANdt  savedBFlist saveAN_spikesOrProbability saveMAPparamsName...
rmeddis@38:     savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ...
rmeddis@38:     DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
rmeddis@38:     IHCoutput ANprobRateOutput ANoutput savePavailable ANtauCas  ...
rmeddis@38:     CNtauGk CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates ...
rmeddis@38:     MOCattenuation
rmeddis@38: 
rmeddis@38: % Normally only ICoutput(logical spike matrix) or ANprobRateOutput will be
rmeddis@38: % needed by the user; so the following will suffice
rmeddis@38: %   global ANdt ICoutput ANprobRateOutput
rmeddis@38: 
rmeddis@38: % Note that sampleRate has not changed from the original function call and
rmeddis@38: %  ANprobRateOutput is sampled at this rate
rmeddis@38: % However ANoutput, CNoutput and IC output are stored as logical
rmeddis@38: %  'spike' matrices using a lower sample rate (see ANdt).
rmeddis@38: 
rmeddis@38: % When AN_spikesOrProbability is set to probability,
rmeddis@38: %  no spike matrices are computed.
rmeddis@38: % When AN_spikesOrProbability is set to 'spikes',
rmeddis@38: %  no probability output is computed
rmeddis@38: 
rmeddis@38: % Efferent control variables are ARattenuation and MOCattenuation
rmeddis@38: %  These are scalars between 1 (no attenuation) and 0.
rmeddis@38: %  They are represented with dt=1/sampleRate (not ANdt)
rmeddis@38: %  They are computed using either AN probability rate output
rmeddis@38: %   or IC (spikes) output as approrpriate.
rmeddis@38: % AR is computed using across channel activity
rmeddis@38: % MOC is computed on a within-channel basis.
rmeddis@38: 
rmeddis@38: if nargin<1
rmeddis@38:     error(' MAP1_14 is not a script but a function that must be called')
rmeddis@38: end
rmeddis@38: 
rmeddis@38: if nargin<6
rmeddis@38:     paramChanges=[];
rmeddis@38: end
rmeddis@38: % Read parameters from MAPparams<***> file in 'parameterStore' folder
rmeddis@38: % Beware, 'BFlist=-1' is a legitimate argument for MAPparams<>
rmeddis@38: %  It means that the calling program allows MAPparams to specify the list
rmeddis@38: cmd=['method=MAPparams' MAPparamsName ...
rmeddis@38:     '(BFlist, sampleRate, 0, paramChanges);'];
rmeddis@38: eval(cmd);
rmeddis@38: BFlist=DRNLParams.nonlinCFs;
rmeddis@38: 
rmeddis@38: % save as global for later plotting if required
rmeddis@38: savedBFlist=BFlist;
rmeddis@38: saveAN_spikesOrProbability=AN_spikesOrProbability;
rmeddis@38: saveMAPparamsName=MAPparamsName;
rmeddis@38: 
rmeddis@38: dt=1/sampleRate;
rmeddis@38: duration=length(inputSignal)/sampleRate;
rmeddis@38: % segmentDuration is specified in parameter file (must be >efferent delay)
rmeddis@38: segmentDuration=method.segmentDuration;
rmeddis@38: segmentLength=round(segmentDuration/ dt);
rmeddis@38: segmentTime=dt*(1:segmentLength); % used in debugging plots
rmeddis@38: 
rmeddis@38: % all spiking activity is computed using longer epochs
rmeddis@38: ANspeedUpFactor=AN_IHCsynapseParams.ANspeedUpFactor;  % e.g.5 times
rmeddis@38: 
rmeddis@38: % inputSignal must be  row vector
rmeddis@38: [r c]=size(inputSignal);
rmeddis@38: if r>c, inputSignal=inputSignal'; end       % transpose
rmeddis@38: % ignore stereo signals
rmeddis@38: inputSignal=inputSignal(1,:);               % drop any second channel
rmeddis@38: savedInputSignal=inputSignal;
rmeddis@38: 
rmeddis@38: % Segment the signal
rmeddis@38: % The sgment length is given but the signal length must be adjusted to be a
rmeddis@38: % multiple of both the segment length and the reduced segmentlength
rmeddis@38: [nSignalRows signalLength]=size(inputSignal);
rmeddis@38: segmentLength=ceil(segmentLength/ANspeedUpFactor)*ANspeedUpFactor;
rmeddis@38: % Make the signal length a whole multiple of the segment length
rmeddis@38: nSignalSegments=ceil(signalLength/segmentLength);
rmeddis@38: padSize=nSignalSegments*segmentLength-signalLength;
rmeddis@38: pad=zeros(nSignalRows,padSize);
rmeddis@38: inputSignal=[inputSignal pad];
rmeddis@38: [ignore signalLength]=size(inputSignal);
rmeddis@38: 
rmeddis@38: % AN (spikes) is computed at a lower sample rate when spikes required
rmeddis@38: %  so it has a reduced segment length (see 'ANspeeUpFactor' above)
rmeddis@38: % AN CN and IC all use this sample interval
rmeddis@38: ANdt=dt*ANspeedUpFactor;
rmeddis@38: reducedSegmentLength=round(segmentLength/ANspeedUpFactor);
rmeddis@38: reducedSignalLength= round(signalLength/ANspeedUpFactor);
rmeddis@38: 
rmeddis@38: %% Initialise with respect to each stage before computing
rmeddis@38: %  by allocating memory,
rmeddis@38: %  by computing constants
rmeddis@38: %  by establishing easy to read variable names
rmeddis@38: % The computations are made in segments and boundary conditions must
rmeddis@38: %  be established and stored. These are found in variables with
rmeddis@38: %  'boundary' or 'bndry' in the name
rmeddis@38: 
rmeddis@38: %% OME ---
rmeddis@38: % external ear resonances
rmeddis@38: OMEexternalResonanceFilters=OMEParams.externalResonanceFilters;
rmeddis@38: [nOMEExtFilters c]=size(OMEexternalResonanceFilters);
rmeddis@38: % details of external (outer ear) resonances
rmeddis@38: OMEgaindBs=OMEexternalResonanceFilters(:,1);
rmeddis@38: OMEgainScalars=10.^(OMEgaindBs/20);
rmeddis@38: OMEfilterOrder=OMEexternalResonanceFilters(:,2);
rmeddis@38: OMElowerCutOff=OMEexternalResonanceFilters(:,3);
rmeddis@38: OMEupperCutOff=OMEexternalResonanceFilters(:,4);
rmeddis@38: % external resonance coefficients
rmeddis@38: ExtFilter_b=cell(nOMEExtFilters,1);
rmeddis@38: ExtFilter_a=cell(nOMEExtFilters,1);
rmeddis@38: for idx=1:nOMEExtFilters
rmeddis@38:     Nyquist=sampleRate/2;
rmeddis@38:     [b, a] = butter(OMEfilterOrder(idx), ...
rmeddis@38:         [OMElowerCutOff(idx) OMEupperCutOff(idx)]...
rmeddis@38:         /Nyquist);
rmeddis@38:     ExtFilter_b{idx}=b;
rmeddis@38:     ExtFilter_a{idx}=a;
rmeddis@38: end
rmeddis@38: OMEExtFilterBndry=cell(2,1);
rmeddis@38: OMEextEarPressure=zeros(1,signalLength); % pressure at tympanic membrane
rmeddis@38: 
rmeddis@38: % pressure to velocity conversion using smoothing filter (50 Hz cutoff)
rmeddis@38: tau=1/(2*pi*50);
rmeddis@38: a1=dt/tau-1; a0=1;
rmeddis@38: b0=1+ a1;
rmeddis@38: TMdisp_b=b0; TMdisp_a=[a0 a1];
rmeddis@38: % figure(9), freqz(TMdisp_b, TMdisp_a)
rmeddis@38: OME_TMdisplacementBndry=[];
rmeddis@38: 
rmeddis@38: % OME high pass (simulates poor low frequency stapes response)
rmeddis@38: OMEhighPassHighCutOff=OMEParams.OMEstapesLPcutoff;
rmeddis@38: Nyquist=sampleRate/2;
rmeddis@38: [stapesDisp_b,stapesDisp_a] = butter(1, OMEhighPassHighCutOff/Nyquist, 'high');
rmeddis@38: % figure(10), freqz(stapesDisp_b, stapesDisp_a)
rmeddis@38: 
rmeddis@38: OMEhighPassBndry=[];
rmeddis@38: 
rmeddis@38: % OMEampStapes might be reducdant (use OMEParams.stapesScalar)
rmeddis@38: stapesScalar= OMEParams.stapesScalar;
rmeddis@38: 
rmeddis@38: % Acoustic reflex
rmeddis@38: efferentDelayPts=round(OMEParams.ARdelay/dt);
rmeddis@38: % smoothing filter
rmeddis@38: a1=dt/OMEParams.ARtau-1; a0=1;
rmeddis@38: b0=1+ a1;
rmeddis@38: ARfilt_b=b0; ARfilt_a=[a0 a1];
rmeddis@38: 
rmeddis@38: ARattenuation=ones(1,signalLength);
rmeddis@38: ARrateThreshold=OMEParams.ARrateThreshold; % may not be used
rmeddis@38: ARrateToAttenuationFactor=OMEParams.rateToAttenuationFactor;
rmeddis@38: ARrateToAttenuationFactorProb=OMEParams.rateToAttenuationFactorProb;
rmeddis@38: ARboundary=[];
rmeddis@38: ARboundaryProb=0;
rmeddis@38: 
rmeddis@38: % save complete OME record (stapes displacement)
rmeddis@38: OMEoutput=zeros(1,signalLength);
rmeddis@38: TMoutput=zeros(1,signalLength);
rmeddis@38: 
rmeddis@38: %% BM ---
rmeddis@38: % BM is represented as a list of locations identified by BF
rmeddis@38: DRNL_BFs=BFlist;
rmeddis@38: nBFs= length(DRNL_BFs);
rmeddis@38: 
rmeddis@38: % DRNLchannelParameters=DRNLParams.channelParameters;
rmeddis@38: DRNLresponse= zeros(nBFs, segmentLength);
rmeddis@38: 
rmeddis@38: MOCrateToAttenuationFactor=DRNLParams.rateToAttenuationFactor;
rmeddis@38: rateToAttenuationFactorProb=DRNLParams.rateToAttenuationFactorProb;
rmeddis@38: MOCrateThresholdProb=DRNLParams.MOCrateThresholdProb;
rmeddis@38: 
rmeddis@38: % smoothing filter for MOC
rmeddis@38: a1=dt/DRNLParams.MOCtau-1; a0=1;
rmeddis@38: b0=1+ a1;
rmeddis@38: MOCfilt_b=b0; MOCfilt_a=[a0 a1];
rmeddis@38: % figure(9), freqz(stapesDisp_b, stapesDisp_a)
rmeddis@38: MOCboundary=cell(nBFs,1);
rmeddis@38: MOCprobBoundary=cell(nBFs,1);
rmeddis@38: 
rmeddis@38: MOCattSegment=zeros(nBFs,reducedSegmentLength);
rmeddis@38: MOCattenuation=ones(nBFs,signalLength);
rmeddis@38: 
rmeddis@38: % if DRNLParams.a>0
rmeddis@38: %     DRNLcompressionThreshold=10^((1/(1-DRNLParams.c))* ...
rmeddis@38: %     log10(DRNLParams.b/DRNLParams.a));
rmeddis@38: % else
rmeddis@38: %     DRNLcompressionThreshold=inf;
rmeddis@38: % end
rmeddis@38: DRNLcompressionThreshold=DRNLParams.cTh;
rmeddis@38: DRNLlinearOrder= DRNLParams.linOrder;
rmeddis@38: DRNLnonlinearOrder= DRNLParams.nonlinOrder;
rmeddis@38: 
rmeddis@38: DRNLa=DRNLParams.a;
rmeddis@38: DRNLb=DRNLParams.b;
rmeddis@38: DRNLc=DRNLParams.c;
rmeddis@38: linGAIN=DRNLParams.g;
rmeddis@38: %
rmeddis@38: % gammatone filter coefficients for linear pathway
rmeddis@38: bw=DRNLParams.linBWs';
rmeddis@38: phi = 2 * pi * bw * dt;
rmeddis@38: cf=DRNLParams.linCFs';
rmeddis@38: theta = 2 * pi * cf * dt;
rmeddis@38: cos_theta = cos(theta);
rmeddis@38: sin_theta = sin(theta);
rmeddis@38: alpha = -exp(-phi).* cos_theta;
rmeddis@38: b0 = ones(nBFs,1);
rmeddis@38: b1 = 2 * alpha;
rmeddis@38: b2 = exp(-2 * phi);
rmeddis@38: z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i;
rmeddis@38: z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i;
rmeddis@38: z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i;
rmeddis@38: tf = (z2 + z3) ./ z1;
rmeddis@38: a0 = abs(tf);
rmeddis@38: a1 = alpha .* a0;
rmeddis@38: GTlin_a = [b0, b1, b2];
rmeddis@38: GTlin_b = [a0, a1];
rmeddis@38: GTlinOrder=DRNLlinearOrder;
rmeddis@38: GTlinBdry=cell(nBFs,GTlinOrder);
rmeddis@38: 
rmeddis@38: % nonlinear gammatone filter coefficients 
rmeddis@38: bw=DRNLParams.nlBWs';
rmeddis@38: phi = 2 * pi * bw * dt;
rmeddis@38: cf=DRNLParams.nonlinCFs';
rmeddis@38: theta = 2 * pi * cf * dt;
rmeddis@38: cos_theta = cos(theta);
rmeddis@38: sin_theta = sin(theta);
rmeddis@38: alpha = -exp(-phi).* cos_theta;
rmeddis@38: b0 = ones(nBFs,1);
rmeddis@38: b1 = 2 * alpha;
rmeddis@38: b2 = exp(-2 * phi);
rmeddis@38: z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i;
rmeddis@38: z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i;
rmeddis@38: z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i;
rmeddis@38: tf = (z2 + z3) ./ z1;
rmeddis@38: a0 = abs(tf);
rmeddis@38: a1 = alpha .* a0;
rmeddis@38: GTnonlin_a = [b0, b1, b2];
rmeddis@38: GTnonlin_b = [a0, a1];
rmeddis@38: GTnonlinOrder=DRNLnonlinearOrder;
rmeddis@38: GTnonlinBdry1=cell(nBFs, GTnonlinOrder);
rmeddis@38: GTnonlinBdry2=cell(nBFs, GTnonlinOrder);
rmeddis@38: 
rmeddis@38: % complete BM record (BM displacement)
rmeddis@38: DRNLoutput=zeros(nBFs, signalLength);
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% IHC ---
rmeddis@38: % IHC cilia activity and receptor potential
rmeddis@38: % viscous coupling between BM and stereocilia displacement
rmeddis@38: % Nyquist=sampleRate/2;
rmeddis@38: % IHCcutoff=1/(2*pi*IHC_cilia_RPParams.tc);
rmeddis@38: % [IHCciliaFilter_b,IHCciliaFilter_a]=...
rmeddis@38: %     butter(1, IHCcutoff/Nyquist, 'high');
rmeddis@38: a1=dt/IHC_cilia_RPParams.tc-1; a0=1;
rmeddis@38: b0=1+ a1;
rmeddis@38: % high pass (i.e. low pass reversed)
rmeddis@38: IHCciliaFilter_b=[a0 a1]; IHCciliaFilter_a=b0;
rmeddis@38: % figure(9), freqz(IHCciliaFilter_b, IHCciliaFilter_a)
rmeddis@38: 
rmeddis@38: IHCciliaBndry=cell(nBFs,1);
rmeddis@38: 
rmeddis@38: % IHC apical conductance (Boltzman function)
rmeddis@38: IHC_C= IHC_cilia_RPParams.C;
rmeddis@38: IHCu0= IHC_cilia_RPParams.u0;
rmeddis@38: IHCu1= IHC_cilia_RPParams.u1;
rmeddis@38: IHCs0= IHC_cilia_RPParams.s0;
rmeddis@38: IHCs1= IHC_cilia_RPParams.s1;
rmeddis@38: IHCGmax= IHC_cilia_RPParams.Gmax;
rmeddis@38: IHCGa= IHC_cilia_RPParams.Ga; % (leakage)
rmeddis@38: 
rmeddis@38: IHCGu0 = IHCGa+IHCGmax./(1+exp(IHCu0/IHCs0).*(1+exp(IHCu1/IHCs1)));
rmeddis@38: IHCrestingCiliaCond=IHCGu0;
rmeddis@38: 
rmeddis@38: % Receptor potential
rmeddis@38: IHC_Cab= IHC_cilia_RPParams.Cab;
rmeddis@38: IHC_Gk= IHC_cilia_RPParams.Gk;
rmeddis@38: IHC_Et= IHC_cilia_RPParams.Et;
rmeddis@38: IHC_Ek= IHC_cilia_RPParams.Ek;
rmeddis@38: IHC_Ekp= IHC_Ek+IHC_Et*IHC_cilia_RPParams.Rpc;
rmeddis@38: 
rmeddis@38: IHCrestingV= (IHC_Gk*IHC_Ekp+IHCGu0*IHC_Et)/(IHCGu0+IHC_Gk);
rmeddis@38: 
rmeddis@38: IHC_Vnow= IHCrestingV*ones(nBFs,1); % initial voltage
rmeddis@38: IHC_RP= zeros(nBFs,segmentLength);
rmeddis@38: 
rmeddis@38: % complete record of IHC receptor potential (V)
rmeddis@38: IHCciliaDisplacement= zeros(nBFs,segmentLength);
rmeddis@38: IHCoutput= zeros(nBFs,signalLength);
rmeddis@38: IHC_cilia_output= zeros(nBFs,signalLength);
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% pre-synapse ---
rmeddis@38: % Each BF is replicated using a different fiber type to make a 'channel'
rmeddis@38: % The number of channels is nBFs x nANfiberTypes
rmeddis@38: % Fiber types are specified in terms of tauCa
rmeddis@38: nANfiberTypes= length(IHCpreSynapseParams.tauCa);
rmeddis@38: ANtauCas= IHCpreSynapseParams.tauCa;
rmeddis@38: nANchannels= nANfiberTypes*nBFs;
rmeddis@38: synapticCa= zeros(nANchannels,segmentLength);
rmeddis@38: 
rmeddis@38: % Calcium control (more calcium, greater release rate)
rmeddis@38: ECa=IHCpreSynapseParams.ECa;
rmeddis@38: gamma=IHCpreSynapseParams.gamma;
rmeddis@38: beta=IHCpreSynapseParams.beta;
rmeddis@38: tauM=IHCpreSynapseParams.tauM;
rmeddis@38: mICa=zeros(nANchannels,segmentLength);
rmeddis@38: GmaxCa=IHCpreSynapseParams.GmaxCa;
rmeddis@38: synapse_z= IHCpreSynapseParams.z;
rmeddis@38: synapse_power=IHCpreSynapseParams.power;
rmeddis@38: 
rmeddis@38: % tauCa vector is established across channels to allow vectorization
rmeddis@38: %  (one tauCa per channel). Do not confuse with ANtauCas (one pre fiber type)
rmeddis@38: tauCa=repmat(ANtauCas, nBFs,1);
rmeddis@38: tauCa=reshape(tauCa, nANchannels, 1);
rmeddis@38: 
rmeddis@38: % presynapse startup values (vectors, length:nANchannels)
rmeddis@38: % proportion (0 - 1) of Ca channels open at IHCrestingV
rmeddis@38: mICaCurrent=((1+beta^-1 * exp(-gamma*IHCrestingV))^-1)...
rmeddis@38:     *ones(nBFs*nANfiberTypes,1);
rmeddis@38: % corresponding startup currents
rmeddis@38: ICaCurrent= (GmaxCa*mICaCurrent.^3) * (IHCrestingV-ECa);
rmeddis@38: CaCurrent= ICaCurrent.*tauCa;
rmeddis@38: 
rmeddis@38: % vesicle release rate at startup (one per channel)
rmeddis@38: % kt0 is used only at initialisation
rmeddis@38: kt0= -synapse_z * CaCurrent.^synapse_power;
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% AN ---
rmeddis@38: % each row of the AN matrices represents one AN fiber
rmeddis@38: % The results computed either for probabiities *or* for spikes (not both)
rmeddis@38: % Spikes are necessary if CN and IC are to be computed
rmeddis@38: nFibersPerChannel= AN_IHCsynapseParams.numFibers;
rmeddis@38: nANfibers= nANchannels*nFibersPerChannel;
rmeddis@38: AN_refractory_period= AN_IHCsynapseParams.refractory_period;
rmeddis@38: 
rmeddis@38: y=AN_IHCsynapseParams.y;
rmeddis@38: l=AN_IHCsynapseParams.l;
rmeddis@38: x=AN_IHCsynapseParams.x;
rmeddis@38: r=AN_IHCsynapseParams.r;
rmeddis@38: M=round(AN_IHCsynapseParams.M);
rmeddis@38: 
rmeddis@38: % probability            (NB initial 'P' on everything)
rmeddis@38: PAN_ydt = repmat(AN_IHCsynapseParams.y*dt, nANchannels,1);
rmeddis@38: PAN_ldt = repmat(AN_IHCsynapseParams.l*dt, nANchannels,1);
rmeddis@38: PAN_xdt = repmat(AN_IHCsynapseParams.x*dt, nANchannels,1);
rmeddis@38: PAN_rdt = repmat(AN_IHCsynapseParams.r*dt, nANchannels,1);
rmeddis@38: PAN_rdt_plus_ldt = PAN_rdt + PAN_ldt;
rmeddis@38: PAN_M=round(AN_IHCsynapseParams.M);
rmeddis@38: 
rmeddis@38: % compute starting values
rmeddis@38: Pcleft    = kt0* y* M ./ (y*(l+r)+ kt0* l);
rmeddis@38: Pavailable    = Pcleft*(l+r)./kt0;
rmeddis@38: Preprocess    = Pcleft*r/x; % canbe fractional
rmeddis@38: 
rmeddis@38: ANprobability=zeros(nANchannels,segmentLength);
rmeddis@38: ANprobRateOutput=zeros(nANchannels,signalLength);
rmeddis@38: lengthAbsRefractoryP= round(AN_refractory_period/dt);
rmeddis@38: cumANnotFireProb=ones(nANchannels,signalLength);
rmeddis@38: % special variables for monitoring synaptic cleft (specialists only)
rmeddis@38: savePavailableSeg=zeros(nANchannels,segmentLength);
rmeddis@38: savePavailable=zeros(nANchannels,signalLength);
rmeddis@38: 
rmeddis@38: % spikes     % !  !  !    ! !        !   !  !
rmeddis@38: lengthAbsRefractory= round(AN_refractory_period/ANdt);
rmeddis@38: 
rmeddis@38: AN_ydt= repmat(AN_IHCsynapseParams.y*ANdt, nANfibers,1);
rmeddis@38: AN_ldt= repmat(AN_IHCsynapseParams.l*ANdt, nANfibers,1);
rmeddis@38: AN_xdt= repmat(AN_IHCsynapseParams.x*ANdt, nANfibers,1);
rmeddis@38: AN_rdt= repmat(AN_IHCsynapseParams.r*ANdt, nANfibers,1);
rmeddis@38: AN_rdt_plus_ldt= AN_rdt + AN_ldt;
rmeddis@38: AN_M= round(AN_IHCsynapseParams.M);
rmeddis@38: 
rmeddis@38: % kt0  is initial release rate
rmeddis@38: % Establish as a vector (length=channel x number of fibers)
rmeddis@38: kt0= repmat(kt0', nFibersPerChannel, 1);
rmeddis@38: kt0=reshape(kt0, nANfibers,1);
rmeddis@38: 
rmeddis@38: % starting values for reservoirs
rmeddis@38: AN_cleft    = kt0* y* M ./ (y*(l+r)+ kt0* l);
rmeddis@38: AN_available    = round(AN_cleft*(l+r)./kt0); %must be integer
rmeddis@38: AN_reprocess    = AN_cleft*r/x;
rmeddis@38: 
rmeddis@38: % output is in a logical array spikes = 1/ 0.
rmeddis@38: ANspikes= false(nANfibers,reducedSegmentLength);
rmeddis@38: ANoutput= false(nANfibers,reducedSignalLength);
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% CN (first brain stem nucleus - could be any subdivision of CN)
rmeddis@38: % Input to a CN neuorn is a random selection of AN fibers within a channel
rmeddis@38: %  The number of AN fibers used is ANfibersFanInToCN
rmeddis@38: % CNtauGk (Potassium time constant) determines the rate of firing of
rmeddis@38: %  the unit when driven hard by a DC input (not normally >350 sp/s)
rmeddis@38: % If there is more than one value, everything is replicated accordingly
rmeddis@38: 
rmeddis@38: ANavailableFibersPerChan=AN_IHCsynapseParams.numFibers;
rmeddis@38: ANfibersFanInToCN=MacGregorMultiParams.fibersPerNeuron;
rmeddis@38: 
rmeddis@38: CNtauGk=MacGregorMultiParams.tauGk; % row vector of CN types (by tauGk)
rmeddis@38: nCNtauGk=length(CNtauGk);
rmeddis@38: 
rmeddis@38: % the total number of 'channels' is now greater
rmeddis@38: % 'channel' is defined as collections of units with the same parameters
rmeddis@38: %  i.e. same BF, same ANtau, same CNtauGk
rmeddis@38: nCNchannels=nANchannels*nCNtauGk;
rmeddis@38: 
rmeddis@38: nCNneuronsPerChannel=MacGregorMultiParams.nNeuronsPerBF;
rmeddis@38: tauGk=repmat(CNtauGk, nCNneuronsPerChannel,1);
rmeddis@38: tauGk=reshape(tauGk,nCNneuronsPerChannel*nCNtauGk,1);
rmeddis@38: 
rmeddis@38: % Now the number of neurons has been increased
rmeddis@38: nCNneurons=nCNneuronsPerChannel*nCNchannels;
rmeddis@38: CNmembranePotential=zeros(nCNneurons,reducedSegmentLength);
rmeddis@38: 
rmeddis@38: % establish which ANfibers (by name) feed into which CN nuerons
rmeddis@38: CNinputfiberLists=zeros(nANchannels*nCNneuronsPerChannel, ANfibersFanInToCN);
rmeddis@38: unitNo=1;
rmeddis@38: for ch=1:nANchannels
rmeddis@38:     % Each channel contains a number of units =length(listOfFanInValues)
rmeddis@38:     for idx=1:nCNneuronsPerChannel
rmeddis@38:         for idx2=1:nCNtauGk
rmeddis@38:             fibersUsed=(ch-1)*ANavailableFibersPerChan + ...
rmeddis@38:                 ceil(rand(1,ANfibersFanInToCN)* ANavailableFibersPerChan);
rmeddis@38:             CNinputfiberLists(unitNo,:)=fibersUsed;
rmeddis@38:             unitNo=unitNo+1;
rmeddis@38:         end
rmeddis@38:     end
rmeddis@38: end
rmeddis@38: 
rmeddis@38: % input to CN units
rmeddis@38: AN_PSTH=zeros(nCNneurons,reducedSegmentLength);
rmeddis@38: 
rmeddis@38: % Generate CNalphaFunction function
rmeddis@38: %  by which spikes are converted to post-synaptic currents
rmeddis@38: CNdendriteLPfreq= MacGregorMultiParams.dendriteLPfreq;
rmeddis@38: CNcurrentPerSpike=MacGregorMultiParams.currentPerSpike;
rmeddis@38: CNspikeToCurrentTau=1/(2*pi*CNdendriteLPfreq);
rmeddis@38: t=ANdt:ANdt:5*CNspikeToCurrentTau;
rmeddis@38: CNalphaFunction= (1 / ...
rmeddis@38:     CNspikeToCurrentTau)*t.*exp(-t /CNspikeToCurrentTau);
rmeddis@38: CNalphaFunction=CNalphaFunction*CNcurrentPerSpike;
rmeddis@38: 
rmeddis@38: % figure(98), plot(t,CNalphaFunction)
rmeddis@38: % working memory for implementing convolution
rmeddis@38: 
rmeddis@38: CNcurrentTemp=...
rmeddis@38:     zeros(nCNneurons,reducedSegmentLength+length(CNalphaFunction)-1);
rmeddis@38: % trailing alphas are parts of humps carried forward to the next segment
rmeddis@38: CNtrailingAlphas=zeros(nCNneurons,length(CNalphaFunction));
rmeddis@38: 
rmeddis@38: CN_tauM=MacGregorMultiParams.tauM;
rmeddis@38: CN_tauTh=MacGregorMultiParams.tauTh;
rmeddis@38: CN_cap=MacGregorMultiParams.Cap;
rmeddis@38: CN_c=MacGregorMultiParams.c;
rmeddis@38: CN_b=MacGregorMultiParams.dGkSpike;
rmeddis@38: CN_Ek=MacGregorMultiParams.Ek;
rmeddis@38: CN_Eb= MacGregorMultiParams.Eb;
rmeddis@38: CN_Er=MacGregorMultiParams.Er;
rmeddis@38: CN_Th0= MacGregorMultiParams.Th0;
rmeddis@38: CN_E= zeros(nCNneurons,1);
rmeddis@38: CN_Gk= zeros(nCNneurons,1);
rmeddis@38: CN_Th= MacGregorMultiParams.Th0*ones(nCNneurons,1);
rmeddis@38: CN_Eb=CN_Eb.*ones(nCNneurons,1);
rmeddis@38: CN_Er=CN_Er.*ones(nCNneurons,1);
rmeddis@38: CNtimeSinceLastSpike=zeros(nCNneurons,1);
rmeddis@38: % tauGk is the main distinction between neurons
rmeddis@38: %  in fact they are all the same in the standard model
rmeddis@38: tauGk=repmat(tauGk,nANchannels,1);
rmeddis@38: 
rmeddis@38: CNoutput=false(nCNneurons,reducedSignalLength);
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% MacGregor (IC - second nucleus) --------
rmeddis@38: nICcells=nANchannels*nCNtauGk;  % one cell per channel
rmeddis@38: CN_PSTH=zeros(nICcells ,reducedSegmentLength);
rmeddis@38: 
rmeddis@38: % ICspikeWidth=0.00015;   % this may need revisiting
rmeddis@38: % epochsPerSpike=round(ICspikeWidth/ANdt);
rmeddis@38: % if epochsPerSpike<1
rmeddis@38: %     error(['MacGregorMulti: sample rate too low to support ' ...
rmeddis@38: %         num2str(ICspikeWidth*1e6) '  microsec spikes']);
rmeddis@38: % end
rmeddis@38: 
rmeddis@38: % short names
rmeddis@38: IC_tauM=MacGregorParams.tauM;
rmeddis@38: IC_tauGk=MacGregorParams.tauGk;
rmeddis@38: IC_tauTh=MacGregorParams.tauTh;
rmeddis@38: IC_cap=MacGregorParams.Cap;
rmeddis@38: IC_c=MacGregorParams.c;
rmeddis@38: IC_b=MacGregorParams.dGkSpike;
rmeddis@38: IC_Th0=MacGregorParams.Th0;
rmeddis@38: IC_Ek=MacGregorParams.Ek;
rmeddis@38: IC_Eb= MacGregorParams.Eb;
rmeddis@38: IC_Er=MacGregorParams.Er;
rmeddis@38: 
rmeddis@38: IC_E=zeros(nICcells,1);
rmeddis@38: IC_Gk=zeros(nICcells,1);
rmeddis@38: IC_Th=IC_Th0*ones(nICcells,1);
rmeddis@38: 
rmeddis@38: % Dendritic filtering, all spikes are replaced by CNalphaFunction functions
rmeddis@38: ICdendriteLPfreq= MacGregorParams.dendriteLPfreq;
rmeddis@38: ICcurrentPerSpike=MacGregorParams.currentPerSpike;
rmeddis@38: ICspikeToCurrentTau=1/(2*pi*ICdendriteLPfreq);
rmeddis@38: t=ANdt:ANdt:3*ICspikeToCurrentTau;
rmeddis@38: IC_CNalphaFunction= (ICcurrentPerSpike / ...
rmeddis@38:     ICspikeToCurrentTau)*t.*exp(-t / ICspikeToCurrentTau);
rmeddis@38: % figure(98), plot(t,IC_CNalphaFunction)
rmeddis@38: 
rmeddis@38: % working space for implementing alpha function
rmeddis@38: ICcurrentTemp=...
rmeddis@38:     zeros(nICcells,reducedSegmentLength+length(IC_CNalphaFunction)-1);
rmeddis@38: ICtrailingAlphas=zeros(nICcells, length(IC_CNalphaFunction));
rmeddis@38: 
rmeddis@38: ICfiberTypeRates=zeros(nANfiberTypes,reducedSignalLength);
rmeddis@38: ICoutput=false(nICcells,reducedSignalLength);
rmeddis@38: 
rmeddis@38: ICmembranePotential=zeros(nICcells,reducedSegmentLength);
rmeddis@38: ICmembraneOutput=zeros(nICcells,signalLength);
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% Main program %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%
rmeddis@38: 
rmeddis@38: %  Compute the entire model for each segment
rmeddis@38: segmentStartPTR=1;
rmeddis@38: reducedSegmentPTR=1; % when sampling rate is reduced
rmeddis@38: while segmentStartPTR<signalLength
rmeddis@38:     segmentEndPTR=segmentStartPTR+segmentLength-1;
rmeddis@38:     % shorter segments after speed up.
rmeddis@38:     shorterSegmentEndPTR=reducedSegmentPTR+reducedSegmentLength-1;
rmeddis@38: 
rmeddis@38:     inputPressureSegment=inputSignal...
rmeddis@38:         (:,segmentStartPTR:segmentStartPTR+segmentLength-1);
rmeddis@38: 
rmeddis@38:     % segment debugging plots
rmeddis@38:     % figure(98)
rmeddis@38:     % plot(segmentTime,inputPressureSegment), title('signalSegment')
rmeddis@38: 
rmeddis@38: 
rmeddis@38:     % OME ----------------------
rmeddis@38: 
rmeddis@38:     % OME Stage 1: external resonances. Add to inputSignal pressure wave
rmeddis@38:     y=inputPressureSegment;
rmeddis@38:     for n=1:nOMEExtFilters
rmeddis@38:         % any number of resonances can be used
rmeddis@38:         [x  OMEExtFilterBndry{n}] = ...
rmeddis@38:             filter(ExtFilter_b{n},ExtFilter_a{n},...
rmeddis@38:             inputPressureSegment, OMEExtFilterBndry{n});
rmeddis@38:         x= x* OMEgainScalars(n);
rmeddis@38:         % This is a parallel resonance so add it
rmeddis@38:         y=y+x;
rmeddis@38:     end
rmeddis@38:     inputPressureSegment=y;
rmeddis@38:     OMEextEarPressure(segmentStartPTR:segmentEndPTR)= inputPressureSegment;
rmeddis@38:     
rmeddis@38:     % OME stage 2: convert input pressure (velocity) to
rmeddis@38:     %  tympanic membrane(TM) displacement using low pass filter
rmeddis@38:     [TMdisplacementSegment  OME_TMdisplacementBndry] = ...
rmeddis@38:         filter(TMdisp_b,TMdisp_a,inputPressureSegment, ...
rmeddis@38:         OME_TMdisplacementBndry);
rmeddis@38:     % and save it
rmeddis@38:     TMoutput(segmentStartPTR:segmentEndPTR)= TMdisplacementSegment;
rmeddis@38: 
rmeddis@38:     % OME stage 3: middle ear high pass effect to simulate stapes inertia
rmeddis@38:     [stapesDisplacement  OMEhighPassBndry] = ...
rmeddis@38:         filter(stapesDisp_b,stapesDisp_a,TMdisplacementSegment, ...
rmeddis@38:         OMEhighPassBndry);
rmeddis@38: 
rmeddis@38:     % OME stage 4:  apply stapes scalar
rmeddis@38:     stapesDisplacement=stapesDisplacement*stapesScalar;
rmeddis@38: 
rmeddis@38:     % OME stage 5:    acoustic reflex stapes attenuation
rmeddis@38:     %  Attenuate the TM response using feedback from LSR fiber activity
rmeddis@38:     if segmentStartPTR>efferentDelayPts
rmeddis@38:         stapesDisplacement= stapesDisplacement.*...
rmeddis@38:             ARattenuation(segmentStartPTR-efferentDelayPts:...
rmeddis@38:             segmentEndPTR-efferentDelayPts);
rmeddis@38:     end
rmeddis@38: 
rmeddis@38:     % segment debugging plots
rmeddis@38:     % figure(98)
rmeddis@38:     % plot(segmentTime, stapesDisplacement), title ('stapesDisplacement')
rmeddis@38: 
rmeddis@38:     % and save
rmeddis@38:     OMEoutput(segmentStartPTR:segmentEndPTR)= stapesDisplacement;
rmeddis@38: 
rmeddis@38: 
rmeddis@38:     %% BM ------------------------------
rmeddis@38:     % Each location is computed separately
rmeddis@38:     for BFno=1:nBFs
rmeddis@38: 
rmeddis@38:         %            *linear* path
rmeddis@38:         linOutput = stapesDisplacement * linGAIN;  % linear gain
rmeddis@38:        
rmeddis@38:         for order = 1 : GTlinOrder
rmeddis@38:             [linOutput GTlinBdry{BFno,order}] = ...
rmeddis@38:                 filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, ...
rmeddis@38:                 GTlinBdry{BFno,order});
rmeddis@38:         end
rmeddis@38: 
rmeddis@38:         %           *nonLinear* path
rmeddis@38:         % efferent attenuation (0 <> 1)
rmeddis@38:         if segmentStartPTR>efferentDelayPts
rmeddis@38:             MOC=MOCattenuation(BFno, segmentStartPTR-efferentDelayPts:...
rmeddis@38:                 segmentEndPTR-efferentDelayPts);
rmeddis@38:         else    % no MOC available yet
rmeddis@38:             MOC=ones(1, segmentLength);
rmeddis@38:         end
rmeddis@38:         % apply MOC to nonlinear input function       
rmeddis@38:         nonlinOutput=stapesDisplacement.* MOC;
rmeddis@38: 
rmeddis@38:         %       first gammatone filter (nonlin path)
rmeddis@38:         for order = 1 : GTnonlinOrder
rmeddis@38:             [nonlinOutput GTnonlinBdry1{BFno,order}] = ...
rmeddis@38:                 filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ...
rmeddis@38:                 nonlinOutput, GTnonlinBdry1{BFno,order});
rmeddis@38:         end
rmeddis@38:         
rmeddis@38: %         %    original   broken stick instantaneous compression
rmeddis@38: %         y= nonlinOutput.* DRNLa;  % linear section.
rmeddis@38: %         % compress parts of the signal above the compression threshold
rmeddis@38: %         abs_x = abs(nonlinOutput);
rmeddis@38: %         idx=find(abs_x>DRNLcompressionThreshold);
rmeddis@38: %         if ~isempty(idx)>0
rmeddis@38: %             y(idx)=sign(y(idx)).* (DRNLb*abs_x(idx).^DRNLc);
rmeddis@38: %         end
rmeddis@38: %         nonlinOutput=y;
rmeddis@38: 
rmeddis@38:         
rmeddis@38:         %   new broken stick instantaneous compression
rmeddis@38:         y= nonlinOutput.* DRNLa;  % linear section attenuation/gain.
rmeddis@38:         % compress parts of the signal above the compression threshold
rmeddis@38: %         holdY=y;
rmeddis@38:         abs_y = abs(y);
rmeddis@38:         idx=find(abs_y>DRNLcompressionThreshold);
rmeddis@38:         if ~isempty(idx)>0
rmeddis@38: %             y(idx)=sign(y(idx)).* (DRNLcompressionThreshold + ...
rmeddis@38: %                 (abs_y(idx)-DRNLcompressionThreshold).^DRNLc);
rmeddis@38:             y(idx)=sign(y(idx)).* (DRNLcompressionThreshold + ...
rmeddis@38:                 (abs_y(idx)-DRNLcompressionThreshold)*DRNLc);
rmeddis@38:         end
rmeddis@38:         nonlinOutput=y;
rmeddis@38:        
rmeddis@38: % % Boltzmann compression
rmeddis@38: % y=(nonlinOutput * DRNLa*100000);
rmeddis@38: % holdY=y;
rmeddis@38: % y=abs(y);
rmeddis@38: % s=10; u=0.0;
rmeddis@38: % x=1./(1+exp(-(y-u)/s))-0.5;
rmeddis@38: % nonlinOutput=sign(nonlinOutput).*x/10000;
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %         if segmentStartPTR==10*segmentLength+1
rmeddis@38: %             figure(90)
rmeddis@38: %         plot(holdY,'b'), hold on
rmeddis@38: %         plot(nonlinOutput, 'r'), hold off
rmeddis@38: %         ylim([-1e-5 1e-5])
rmeddis@38: %         pause(1)
rmeddis@38: %         end
rmeddis@38: 
rmeddis@38:         %       second filter removes distortion products
rmeddis@38:         for order = 1 : GTnonlinOrder
rmeddis@38:             [ nonlinOutput GTnonlinBdry2{BFno,order}] = ...
rmeddis@38:                 filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ...
rmeddis@38:                 nonlinOutput, GTnonlinBdry2{BFno,order});
rmeddis@38:         end
rmeddis@38: 
rmeddis@38:         %  combine the two paths to give the DRNL displacement
rmeddis@38:         DRNLresponse(BFno,:)=linOutput+nonlinOutput;
rmeddis@38: %         disp(num2str(max(linOutput)))
rmeddis@38:     end % BF
rmeddis@38: 
rmeddis@38:     % segment debugging plots
rmeddis@38:     % figure(98)
rmeddis@38:     %     if size(DRNLresponse,1)>3
rmeddis@38:     %         imagesc(DRNLresponse)  % matrix display
rmeddis@38:     %         title('DRNLresponse'); % single or double channel response
rmeddis@38:     %     else
rmeddis@38:     %         plot(segmentTime, DRNLresponse)
rmeddis@38:     %     end
rmeddis@38: 
rmeddis@38:     % and save it
rmeddis@38:     DRNLoutput(:, segmentStartPTR:segmentEndPTR)= DRNLresponse;
rmeddis@38: 
rmeddis@38: 
rmeddis@38:     %% IHC ------------------------------------
rmeddis@38:     %  BM displacement to IHCciliaDisplacement is  a high-pass filter
rmeddis@38:     %   because of viscous coupling
rmeddis@38:     for idx=1:nBFs
rmeddis@38:         [IHCciliaDisplacement(idx,:)  IHCciliaBndry{idx}] = ...
rmeddis@38:             filter(IHCciliaFilter_b,IHCciliaFilter_a, ...
rmeddis@38:             DRNLresponse(idx,:), IHCciliaBndry{idx});
rmeddis@38:     end
rmeddis@38:     
rmeddis@38:     % apply scalar
rmeddis@38:     IHCciliaDisplacement=IHCciliaDisplacement* IHC_C;
rmeddis@38: 
rmeddis@38:     % and save it
rmeddis@38:     IHC_cilia_output(:,segmentStartPTR:segmentStartPTR+segmentLength-1)=...
rmeddis@38:         IHCciliaDisplacement;
rmeddis@38: 
rmeddis@38:     % compute apical conductance
rmeddis@38:     G=IHCGmax./(1+exp(-(IHCciliaDisplacement-IHCu0)/IHCs0).*...
rmeddis@38:         (1+exp(-(IHCciliaDisplacement-IHCu1)/IHCs1)));
rmeddis@38:     Gu=G + IHCGa;
rmeddis@38: 
rmeddis@38:     % Compute receptor potential
rmeddis@38:     for idx=1:segmentLength
rmeddis@38:         IHC_Vnow=IHC_Vnow+ (-Gu(:, idx).*(IHC_Vnow-IHC_Et)-...
rmeddis@38:             IHC_Gk*(IHC_Vnow-IHC_Ekp))*  dt/IHC_Cab;
rmeddis@38:         IHC_RP(:,idx)=IHC_Vnow;
rmeddis@38:     end
rmeddis@38: 
rmeddis@38:     % segment debugging plots
rmeddis@38:     %     if size(IHC_RP,1)>3
rmeddis@38:     %         surf(IHC_RP), shading interp, title('IHC_RP')
rmeddis@38:     %     else
rmeddis@38:     %         plot(segmentTime, IHC_RP)
rmeddis@38:     %     end
rmeddis@38: 
rmeddis@38:     % and save it
rmeddis@38:     IHCoutput(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=IHC_RP;
rmeddis@38: 
rmeddis@38: 
rmeddis@38:     %% synapse -----------------------------
rmeddis@38:     % Compute the vesicle release rate for each fiber type at each BF
rmeddis@38:     % replicate IHC_RP for each fiber type
rmeddis@38:     Vsynapse=repmat(IHC_RP, nANfiberTypes,1);
rmeddis@38: 
rmeddis@38:     % look-up table of target fraction channels open for a given IHC_RP
rmeddis@38:     mICaINF=    1./( 1 + exp(-gamma  * Vsynapse)  /beta);
rmeddis@38:     % fraction of channel open - apply time constant
rmeddis@38:     for idx=1:segmentLength
rmeddis@38:         % mICaINF is the current 'target' value of mICa
rmeddis@38:         mICaCurrent=mICaCurrent+(mICaINF(:,idx)-mICaCurrent)*dt./tauM;
rmeddis@38:         mICa(:,idx)=mICaCurrent;
rmeddis@38:     end
rmeddis@38: 
rmeddis@38:     ICa=   (GmaxCa* mICa.^3) .* (Vsynapse- ECa);
rmeddis@38: 
rmeddis@38:     for idx=1:segmentLength
rmeddis@38:         CaCurrent=CaCurrent +  ICa(:,idx)*dt - CaCurrent*dt./tauCa;
rmeddis@38:         synapticCa(:,idx)=CaCurrent;
rmeddis@38:     end
rmeddis@38:     synapticCa=-synapticCa; % treat synapticCa as positive substance
rmeddis@38: 
rmeddis@38:     % NB vesicleReleaseRate is /s and is independent of dt
rmeddis@38:     vesicleReleaseRate = synapse_z * synapticCa.^synapse_power; % rate
rmeddis@38: 
rmeddis@38:     % segment debugging plots
rmeddis@38:     %     if size(vesicleReleaseRate,1)>3
rmeddis@38:     %         surf(vesicleReleaseRate), shading interp, title('vesicleReleaseRate')
rmeddis@38:     %     else
rmeddis@38:     %         plot(segmentTime, vesicleReleaseRate)
rmeddis@38:     %     end
rmeddis@38: 
rmeddis@38: 
rmeddis@38:     %% AN
rmeddis@38:     switch AN_spikesOrProbability
rmeddis@38:         case 'probability'
rmeddis@38:             % No refractory effect is applied
rmeddis@38:             for t = 1:segmentLength;
rmeddis@38:                 M_Pq=PAN_M-Pavailable;
rmeddis@38:                 M_Pq(M_Pq<0)=0;
rmeddis@38:                 Preplenish = M_Pq .* PAN_ydt;
rmeddis@38:                 Pejected = Pavailable.* vesicleReleaseRate(:,t)*dt;
rmeddis@38:                 Preprocessed = M_Pq.*Preprocess.* PAN_xdt;
rmeddis@38: 
rmeddis@38:                 ANprobability(:,t)= min(Pejected,1);
rmeddis@38:                 reuptakeandlost= PAN_rdt_plus_ldt .* Pcleft;
rmeddis@38:                 reuptake= PAN_rdt.* Pcleft;
rmeddis@38: 
rmeddis@38:                 Pavailable= Pavailable+ Preplenish- Pejected+ Preprocessed;
rmeddis@38:                 Pcleft= Pcleft + Pejected - reuptakeandlost;
rmeddis@38:                 Preprocess= Preprocess + reuptake - Preprocessed;
rmeddis@38:                 Pavailable(Pavailable<0)=0;
rmeddis@38:                 savePavailableSeg(:,t)=Pavailable;    % synapse tracking
rmeddis@38:                 
rmeddis@38:             end
rmeddis@38: 
rmeddis@38:             % and save it as *rate*
rmeddis@38:             ANrate=ANprobability/dt;
rmeddis@38:             ANprobRateOutput(:, segmentStartPTR:...
rmeddis@38:                 segmentStartPTR+segmentLength-1)=  ANrate;
rmeddis@38:             % monitor synapse contents (only sometimes used)
rmeddis@38:             savePavailable(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=...
rmeddis@38:                 savePavailableSeg;
rmeddis@38:             
rmeddis@38:             %% Apply refractory effect
rmeddis@38:                % the probability of a spike's occurring in the preceding
rmeddis@38:                %  refractory window (t= tnow-refractory period to tnow-)
rmeddis@38:                %    pFired= 1 - II(1-p(t)),
rmeddis@38:                % we need a running account of cumProb=II(1-p(t))
rmeddis@38:                %   cumProb(t)= cumProb(t-1)*(1-p(t))/(1-p(t-refracPeriod))
rmeddis@38:                %   cumProb(0)=0
rmeddis@38:                %   pFired(t)= 1-cumProb(t)
rmeddis@38:                % This gives the fraction of firing events that must be 
rmeddis@38:                %  discounted because of a firing event in the refractory
rmeddis@38:                %  period
rmeddis@38:                %   p(t)= ANprobOutput(t) * pFired(t)
rmeddis@38:                % where ANprobOutput is the uncorrected firing probability
rmeddis@38:                %  based on vesicle release rate
rmeddis@38:                % NB this covers only the absoute refractory period
rmeddis@38:                % not the relative refractory period. To approximate this it
rmeddis@38:                % is necessary to extend the refractory period by 50%
rmeddis@38: 
rmeddis@38: 
rmeddis@38:                 for t = segmentStartPTR:segmentEndPTR;
rmeddis@38:                     if t>1
rmeddis@38:                     ANprobRateOutput(:,t)= ANprobRateOutput(:,t)...
rmeddis@38:                         .* cumANnotFireProb(:,t-1);
rmeddis@38:                     end
rmeddis@38:                     % add recent and remove distant probabilities
rmeddis@38:                     refrac=round(lengthAbsRefractoryP * 1.5);
rmeddis@38:                     if t>refrac
rmeddis@38:                         cumANnotFireProb(:,t)= cumANnotFireProb(:,t-1)...
rmeddis@38:                             .*(1-ANprobRateOutput(:,t)*dt)...
rmeddis@38:                             ./(1-ANprobRateOutput(:,t-refrac)*dt);
rmeddis@38:                     end
rmeddis@38:                 end
rmeddis@38: %                 figure(88), plot(cumANnotFireProb'), title('cumNotFire')
rmeddis@38: %                 figure(89), plot(ANprobRateOutput'), title('ANprobRateOutput')
rmeddis@38: 
rmeddis@38:             %% Estimate efferent effects. ARattenuation (0 <> 1)
rmeddis@38:             %  acoustic reflex
rmeddis@38:             [r c]=size(ANrate);
rmeddis@38:             if r>nBFs % Only if LSR fibers are computed
rmeddis@38:                 ARAttSeg=mean(ANrate(1:nBFs,:),1); %LSR channels are 1:nBF
rmeddis@38:                 % smooth
rmeddis@38:                 [ARAttSeg, ARboundaryProb] = ...
rmeddis@38:                     filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundaryProb);
rmeddis@38:                 ARAttSeg=ARAttSeg-ARrateThreshold;
rmeddis@38:                 ARAttSeg(ARAttSeg<0)=0;   % prevent negative strengths
rmeddis@38:                 ARattenuation(segmentStartPTR:segmentEndPTR)=...
rmeddis@38:                     (1-ARrateToAttenuationFactorProb.* ARAttSeg);
rmeddis@38:             end
rmeddis@38:             %             plot(ARattenuation)
rmeddis@38: 
rmeddis@38:             % MOC attenuation based on within-channel HSR fiber activity
rmeddis@38:             HSRbegins=nBFs*(nANfiberTypes-1)+1;
rmeddis@38:             rates=ANrate(HSRbegins:end,:);
rmeddis@38:             if rateToAttenuationFactorProb<0
rmeddis@38:                 % negative factor implies a fixed attenuation
rmeddis@38:                 MOCattenuation(:,segmentStartPTR:segmentEndPTR)= ...
rmeddis@38:                     ones(size(rates))* -rateToAttenuationFactorProb;
rmeddis@38:             else
rmeddis@38:                 for idx=1:nBFs
rmeddis@38:                     [smoothedRates, MOCprobBoundary{idx}] = ...
rmeddis@38:                         filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ...
rmeddis@38:                         MOCprobBoundary{idx});
rmeddis@38:                     smoothedRates=smoothedRates-MOCrateThresholdProb;
rmeddis@38:                     smoothedRates=max(smoothedRates, 0);
rmeddis@38:                     
rmeddis@38:                     x=(1- smoothedRates* rateToAttenuationFactorProb);
rmeddis@38:                     x=max(x, 10^(-30/20));
rmeddis@38:                     MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= x;
rmeddis@38:                 end
rmeddis@38:             end
rmeddis@38:             MOCattenuation(MOCattenuation<0)=0.001;
rmeddis@38: 
rmeddis@38:             %             plot(MOCattenuation)
rmeddis@38: 
rmeddis@38: 
rmeddis@38:         case 'spikes'
rmeddis@38:             ANtimeCount=0;
rmeddis@38:             % implement speed upt
rmeddis@38:             for t = ANspeedUpFactor:ANspeedUpFactor:segmentLength;
rmeddis@38:                 ANtimeCount=ANtimeCount+1;
rmeddis@38:                 % convert release rate to probabilities
rmeddis@38:                 releaseProb=vesicleReleaseRate(:,t)*ANdt;
rmeddis@38:                 % releaseProb is the release probability per channel
rmeddis@38:                 %  but each channel has many synapses
rmeddis@38:                 releaseProb=repmat(releaseProb',nFibersPerChannel,1);
rmeddis@38:                 releaseProb=reshape(releaseProb, nFibersPerChannel*nANchannels,1);
rmeddis@38: 
rmeddis@38:                 % AN_available=round(AN_available); % vesicles must be integer, (?needed)
rmeddis@38:                 M_q=AN_M- AN_available;     % number of missing vesicles
rmeddis@38:                 M_q(M_q<0)= 0;              % cannot be less than 0
rmeddis@38: 
rmeddis@38:                 % AN_N1 converts probability to discrete events
rmeddis@38:                 %   it considers each event that might occur
rmeddis@38:                 %   (how many vesicles might be released)
rmeddis@38:                 %   and returns a count of how many were released
rmeddis@38: 
rmeddis@38:                 % slow line
rmeddis@38: %                 probabilities= 1-(1-releaseProb).^AN_available;
rmeddis@38:                 probabilities= 1-intpow((1-releaseProb), AN_available);
rmeddis@38:                 ejected= probabilities> rand(length(AN_available),1);
rmeddis@38: 
rmeddis@38:                 reuptakeandlost = AN_rdt_plus_ldt .* AN_cleft;
rmeddis@38:                 reuptake = AN_rdt.* AN_cleft;
rmeddis@38: 
rmeddis@38:                 % slow line
rmeddis@38: %                 probabilities= 1-(1-AN_reprocess.*AN_xdt).^M_q;
rmeddis@38:                 probabilities= 1-intpow((1-AN_reprocess.*AN_xdt), M_q);
rmeddis@38:                 reprocessed= probabilities>rand(length(M_q),1);
rmeddis@38: 
rmeddis@38:                 % slow line
rmeddis@38: %                 probabilities= 1-(1-AN_ydt).^M_q;
rmeddis@38:                  probabilities= 1-intpow((1-AN_ydt), M_q);
rmeddis@38: 
rmeddis@38:                 replenish= probabilities>rand(length(M_q),1);
rmeddis@38: 
rmeddis@38:                 AN_available = AN_available + replenish - ejected ...
rmeddis@38:                     + reprocessed;
rmeddis@38:                 AN_cleft = AN_cleft + ejected - reuptakeandlost;
rmeddis@38:                 AN_reprocess = AN_reprocess + reuptake - reprocessed;
rmeddis@38: 
rmeddis@38:                 % ANspikes is logical record of vesicle release events>0
rmeddis@38:                 ANspikes(:, ANtimeCount)= ejected;
rmeddis@38:             end % t
rmeddis@38: 
rmeddis@38:             % zero any events that are preceded by release events ...
rmeddis@38:             %  within the refractory period
rmeddis@38:             % The refractory period consist of two periods
rmeddis@38:             %  1) the absolute period where no spikes occur
rmeddis@38:             %  2) a relative period where a spike may occur. This relative
rmeddis@38:             %     period is realised as a variable length interval
rmeddis@38:             %     where the length is chosen at random
rmeddis@38:             %     (esentially a linear ramp up)
rmeddis@38: 
rmeddis@38:             % Andreas has a fix for this
rmeddis@38:             for t = 1:ANtimeCount-2*lengthAbsRefractory;
rmeddis@38:                 % identify all spikes across fiber array at time (t)
rmeddis@38:                 % idx is a list of channels where spikes occurred
rmeddis@38:                 % ?? try sparse matrices?
rmeddis@38:                 idx=find(ANspikes(:,t));
rmeddis@38:                 for j=idx  % consider each spike
rmeddis@38:                     % specify variable refractory period
rmeddis@38:                     %  between abs and 2*abs refractory period
rmeddis@38:                     nPointsRefractory=lengthAbsRefractory+...
rmeddis@38:                         round(rand*lengthAbsRefractory);
rmeddis@38:                     % disable spike potential for refractory period
rmeddis@38:                     % set all values in this range to 0
rmeddis@38:                     ANspikes(j,t+1:t+nPointsRefractory)=0;
rmeddis@38:                 end
rmeddis@38:             end  %t
rmeddis@38: 
rmeddis@38:             % segment debugging
rmeddis@38:             % plotInstructions.figureNo=98;
rmeddis@38:             % plotInstructions.displaydt=ANdt;
rmeddis@38:             %  plotInstructions.numPlots=1;
rmeddis@38:             %  plotInstructions.subPlotNo=1;
rmeddis@38:             % UTIL_plotMatrix(ANspikes, plotInstructions);
rmeddis@38: 
rmeddis@38:             % and save it. NB, AN is now on 'speedUp' time
rmeddis@38:             ANoutput(:, reducedSegmentPTR: shorterSegmentEndPTR)=ANspikes;
rmeddis@38: 
rmeddis@38: 
rmeddis@38:             %% CN Macgregor first neucleus -------------------------------
rmeddis@38:             % input is from AN so ANdt is used throughout
rmeddis@38:             % Each CNneuron has a unique set of input fibers selected
rmeddis@38:             %  at random from the available AN fibers (CNinputfiberLists)
rmeddis@38: 
rmeddis@38:             % Create the dendritic current for that neuron
rmeddis@38:             % First get input spikes to this neuron
rmeddis@38:             synapseNo=1;
rmeddis@38:             for ch=1:nCNchannels
rmeddis@38:                 for idx=1:nCNneuronsPerChannel
rmeddis@38:                     % determine candidate fibers for this unit
rmeddis@38:                     fibersUsed=CNinputfiberLists(synapseNo,:);
rmeddis@38:                     % ANpsth has a bin width of ANdt
rmeddis@38:                     %  (just a simple sum across fibers)
rmeddis@38:                     AN_PSTH(synapseNo,:) = ...
rmeddis@38:                         sum(ANspikes(fibersUsed,:), 1);
rmeddis@38:                     synapseNo=synapseNo+1;
rmeddis@38:                 end
rmeddis@38:             end
rmeddis@38: 
rmeddis@38:             % One alpha function per spike
rmeddis@38:             [alphaRows alphaCols]=size(CNtrailingAlphas);
rmeddis@38: 
rmeddis@38:             for unitNo=1:nCNneurons
rmeddis@38:                 CNcurrentTemp(unitNo,:)= ...
rmeddis@38:                     conv2(AN_PSTH(unitNo,:),CNalphaFunction);
rmeddis@38:             end
rmeddis@38: %             disp(['sum(AN_PSTH)= ' num2str(sum(AN_PSTH(1,:)))])
rmeddis@38:             % add post-synaptic current  left over from previous segment
rmeddis@38:             CNcurrentTemp(:,1:alphaCols)=...
rmeddis@38:                 CNcurrentTemp(:,1:alphaCols)+ CNtrailingAlphas;
rmeddis@38: 
rmeddis@38:             % take post-synaptic current for this segment
rmeddis@38:             CNcurrentInput= CNcurrentTemp(:, 1:reducedSegmentLength);
rmeddis@38: %                 disp(['mean(CNcurrentInput)= ' num2str(mean(CNcurrentInput(1,:)))])
rmeddis@38: 
rmeddis@38:             % trailingalphas are the ends of the alpha functions that
rmeddis@38:             % spill over into the next segment
rmeddis@38:             CNtrailingAlphas= ...
rmeddis@38:                 CNcurrentTemp(:, reducedSegmentLength+1:end);
rmeddis@38: 
rmeddis@38:             if CN_c>0
rmeddis@38:                 % variable threshold condition (slow)
rmeddis@38:                 for t=1:reducedSegmentLength
rmeddis@38:                     CNtimeSinceLastSpike=CNtimeSinceLastSpike-ANdt;
rmeddis@38:                     s=CN_E>CN_Th & CNtimeSinceLastSpike<0 ;
rmeddis@38:                     CNtimeSinceLastSpike(s)=0.0005;         % 0.5 ms for sodium spike
rmeddis@38:                     dE =(-CN_E/CN_tauM + ...
rmeddis@38:                         CNcurrentInput(:,t)/CN_cap+(...
rmeddis@38:                         CN_Gk/CN_cap).*(CN_Ek-CN_E))*ANdt;
rmeddis@38:                     dGk=-CN_Gk*ANdt./tauGk + CN_b*s;
rmeddis@38:                     dTh=-(CN_Th-CN_Th0)*ANdt/CN_tauTh + CN_c*s;
rmeddis@38:                     CN_E=CN_E+dE;
rmeddis@38:                     CN_Gk=CN_Gk+dGk;
rmeddis@38:                     CN_Th=CN_Th+dTh;
rmeddis@38:                     CNmembranePotential(:,t)=CN_E+s.*(CN_Eb-CN_E)+CN_Er;
rmeddis@38:                 end
rmeddis@38:             else
rmeddis@38:                 % static threshold (faster)
rmeddis@38:                 E=zeros(1,reducedSegmentLength);
rmeddis@38:                 Gk=zeros(1,reducedSegmentLength);
rmeddis@38:                 ss=zeros(1,reducedSegmentLength);
rmeddis@38:                 for t=1:reducedSegmentLength
rmeddis@38:                     % time of previous spike moves back in time
rmeddis@38:                     CNtimeSinceLastSpike=CNtimeSinceLastSpike-ANdt;
rmeddis@38:                     % action potential if E>threshold
rmeddis@38:                     %  allow time for s to reset between events
rmeddis@38:                     s=CN_E>CN_Th0 & CNtimeSinceLastSpike<0 ;  
rmeddis@38:                     ss(t)=s(1);
rmeddis@38:                     CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike
rmeddis@38:                     dE = (-CN_E/CN_tauM + ...
rmeddis@38:                         CNcurrentInput(:,t)/CN_cap +...
rmeddis@38:                         (CN_Gk/CN_cap).*(CN_Ek-CN_E))*ANdt;
rmeddis@38:                     dGk=-CN_Gk*ANdt./tauGk +CN_b*s;
rmeddis@38:                     CN_E=CN_E+dE;
rmeddis@38:                     CN_Gk=CN_Gk+dGk;
rmeddis@38:                     E(t)=CN_E(1);
rmeddis@38:                     Gk(t)=CN_Gk(1);
rmeddis@38:                     % add spike to CN_E and add resting potential (-60 mV)
rmeddis@38:                     CNmembranePotential(:,t)=CN_E +s.*(CN_Eb-CN_E)+CN_Er;
rmeddis@38:                 end
rmeddis@38:             end
rmeddis@38: %             disp(['CN_E= ' num2str(sum(CN_E(1,:)))])
rmeddis@38: %             disp(['CN_Gk= ' num2str(sum(CN_Gk(1,:)))])
rmeddis@38: %             disp(['CNmembranePotential= ' num2str(sum(CNmembranePotential(1,:)))])
rmeddis@38: %             plot(CNmembranePotential(1,:))
rmeddis@38: 
rmeddis@38: 
rmeddis@38:             % extract spikes.  A spike is a substantial upswing in voltage
rmeddis@38:             CN_spikes=CNmembranePotential> -0.02;
rmeddis@38: %             disp(['CNspikesbefore= ' num2str(sum(sum(CN_spikes)))])
rmeddis@38: 
rmeddis@38:             % now remove any spike that is immediately followed by a spike
rmeddis@38:             % NB 'find' works on columns (whence the transposing)
rmeddis@38:             % for each spike put a zero in the next epoch
rmeddis@38:             CN_spikes=CN_spikes';
rmeddis@38:             idx=find(CN_spikes);
rmeddis@38:             idx=idx(1:end-1);
rmeddis@38:             CN_spikes(idx+1)=0;
rmeddis@38:             CN_spikes=CN_spikes';
rmeddis@38: %             disp(['CNspikes= ' num2str(sum(sum(CN_spikes)))])
rmeddis@38: 
rmeddis@38:             % segment debugging
rmeddis@38:             % plotInstructions.figureNo=98;
rmeddis@38:             % plotInstructions.displaydt=ANdt;
rmeddis@38:             %  plotInstructions.numPlots=1;
rmeddis@38:             %  plotInstructions.subPlotNo=1;
rmeddis@38:             % UTIL_plotMatrix(CN_spikes, plotInstructions);
rmeddis@38: 
rmeddis@38:             % and save it
rmeddis@38:             CNoutput(:, reducedSegmentPTR:shorterSegmentEndPTR)=...
rmeddis@38:                 CN_spikes;
rmeddis@38: 
rmeddis@38: 
rmeddis@38:             %% IC ----------------------------------------------
rmeddis@38:                 %  MacGregor or some other second order neurons
rmeddis@38: 
rmeddis@38:                 % combine CN neurons in same channel (BF x AN tau x CNtau) 
rmeddis@38:                 %  i.e. same BF, same tauCa, same CNtau
rmeddis@38:                 %  to generate inputs to single IC unit
rmeddis@38:                 channelNo=0;
rmeddis@38:                 for idx=1:nCNneuronsPerChannel: ...
rmeddis@38:                         nCNneurons-nCNneuronsPerChannel+1;
rmeddis@38:                     channelNo=channelNo+1;
rmeddis@38:                     CN_PSTH(channelNo,:)=...
rmeddis@38:                         sum(CN_spikes(idx:idx+nCNneuronsPerChannel-1,:));
rmeddis@38:                 end
rmeddis@38: 
rmeddis@38:                 [alphaRows alphaCols]=size(ICtrailingAlphas);
rmeddis@38:                 for ICneuronNo=1:nICcells
rmeddis@38:                     ICcurrentTemp(ICneuronNo,:)= ...
rmeddis@38:                         conv2(CN_PSTH(ICneuronNo,:),  IC_CNalphaFunction);
rmeddis@38:                 end
rmeddis@38: 
rmeddis@38:                 % add the unused current from the previous convolution
rmeddis@38:                 ICcurrentTemp(:,1:alphaCols)=ICcurrentTemp(:,1:alphaCols)...
rmeddis@38:                     + ICtrailingAlphas;
rmeddis@38:                 % take what is required and keep the trailing part for next time
rmeddis@38:                 inputCurrent=ICcurrentTemp(:, 1:reducedSegmentLength);
rmeddis@38:                 ICtrailingAlphas=ICcurrentTemp(:, reducedSegmentLength+1:end);
rmeddis@38: 
rmeddis@38:                 if IC_c==0
rmeddis@38:                     % faster computation when threshold is stable (c==0)
rmeddis@38:                     for t=1:reducedSegmentLength
rmeddis@38:                         s=IC_E>IC_Th0;
rmeddis@38:                         dE = (-IC_E/IC_tauM + inputCurrent(:,t)/IC_cap +...
rmeddis@38:                             (IC_Gk/IC_cap).*(IC_Ek-IC_E))*ANdt;
rmeddis@38:                         dGk=-IC_Gk*ANdt/IC_tauGk +IC_b*s;
rmeddis@38:                         IC_E=IC_E+dE;
rmeddis@38:                         IC_Gk=IC_Gk+dGk;
rmeddis@38:                         ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er;
rmeddis@38:                     end
rmeddis@38:                 else
rmeddis@38:                     %  threshold is changing (IC_c>0; e.g. bushy cell)
rmeddis@38:                     for t=1:reducedSegmentLength
rmeddis@38:                         dE = (-IC_E/IC_tauM + ...
rmeddis@38:                             inputCurrent(:,t)/IC_cap + (IC_Gk/IC_cap)...
rmeddis@38:                             .*(IC_Ek-IC_E))*ANdt;
rmeddis@38:                         IC_E=IC_E+dE;
rmeddis@38:                         s=IC_E>IC_Th;
rmeddis@38:                         ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er;
rmeddis@38:                         dGk=-IC_Gk*ANdt/IC_tauGk +IC_b*s;
rmeddis@38:                         IC_Gk=IC_Gk+dGk;
rmeddis@38: 
rmeddis@38:                         % After a spike, the threshold is raised
rmeddis@38:                         % otherwise it settles to its baseline
rmeddis@38:                         dTh=-(IC_Th-Th0)*ANdt/IC_tauTh +s*IC_c;
rmeddis@38:                         IC_Th=IC_Th+dTh;
rmeddis@38:                     end
rmeddis@38:                 end
rmeddis@38: 
rmeddis@38:                 ICspikes=ICmembranePotential> -0.01;
rmeddis@38:                 % now remove any spike that is immediately followed by a spike
rmeddis@38:                 % NB 'find' works on columns (whence the transposing)
rmeddis@38:                 ICspikes=ICspikes';
rmeddis@38:                 idx=find(ICspikes);
rmeddis@38:                 idx=idx(1:end-1);
rmeddis@38:                 ICspikes(idx+1)=0;
rmeddis@38:                 ICspikes=ICspikes';
rmeddis@38: 
rmeddis@38:                 nCellsPerTau= nICcells/nANfiberTypes;
rmeddis@38:                 firstCell=1;
rmeddis@38:                 lastCell=nCellsPerTau;
rmeddis@38:                 for tauCount=1:nANfiberTypes
rmeddis@38:                     % separate rates according to fiber types
rmeddis@38:                     % currently only the last segment is saved
rmeddis@38:                     ICfiberTypeRates(tauCount, ...
rmeddis@38:                         reducedSegmentPTR:shorterSegmentEndPTR)=...
rmeddis@38:                         sum(ICspikes(firstCell:lastCell, :))...
rmeddis@38:                         /(nCellsPerTau*ANdt);
rmeddis@38:                     firstCell=firstCell+nCellsPerTau;
rmeddis@38:                     lastCell=lastCell+nCellsPerTau;
rmeddis@38:                 end
rmeddis@38:                 
rmeddis@38:                 ICoutput(:,reducedSegmentPTR:shorterSegmentEndPTR)=ICspikes;
rmeddis@38:                 
rmeddis@38:                 % store membrane output on original dt scale
rmeddis@38:                 % do this for single channel models only 
rmeddis@38:                 % and only for the HSR-driven IC cells
rmeddis@38:                 if round(nICcells/length(ANtauCas))==1  % single channel
rmeddis@38:                     % select HSR driven cells
rmeddis@38:                     x= ICmembranePotential(length(ANtauCas),:);
rmeddis@38:                     % restore original dt
rmeddis@38:                     x= repmat(x, ANspeedUpFactor,1);
rmeddis@38:                     x= reshape(x,1,segmentLength);
rmeddis@38:                     if nANfiberTypes>1  % save HSR and LSR
rmeddis@38:                         y=repmat(ICmembranePotential(end,:),...
rmeddis@38:                             ANspeedUpFactor,1);
rmeddis@38:                         y= reshape(y,1,segmentLength);
rmeddis@38:                         x=[x; y];
rmeddis@38:                     end
rmeddis@38:                     ICmembraneOutput(:, segmentStartPTR:segmentEndPTR)= x;
rmeddis@38:                 end
rmeddis@38: 
rmeddis@38:                 % estimate efferent effects.
rmeddis@38:                 % ARis based on LSR units. LSR channels are 1:nBF
rmeddis@38:                 if nANfiberTypes>1  % AR is multi-channel only
rmeddis@38:                     ARAttSeg=sum(ICspikes(1:nBFs,:),1)/ANdt;
rmeddis@38:                     [ARAttSeg, ARboundary] = ...
rmeddis@38:                         filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundary);
rmeddis@38:                     ARAttSeg=ARAttSeg-ARrateThreshold;
rmeddis@38:                     ARAttSeg(ARAttSeg<0)=0;   % prevent negative strengths
rmeddis@38:                     % scale up to dt from ANdt
rmeddis@38:                     x=    repmat(ARAttSeg, ANspeedUpFactor,1);
rmeddis@38:                     x=reshape(x,1,segmentLength);
rmeddis@38:                     ARattenuation(segmentStartPTR:segmentEndPTR)=...
rmeddis@38:                         (1-ARrateToAttenuationFactor* x);
rmeddis@38:                     ARattenuation(ARattenuation<0)=0.001;
rmeddis@38:                 else
rmeddis@38:                     % single channel model; disable AR
rmeddis@38:                     ARattenuation(segmentStartPTR:segmentEndPTR)=...
rmeddis@38:                         ones(1,segmentLength);
rmeddis@38:                 end
rmeddis@38: 
rmeddis@38:                 % MOC attenuation using HSR response only
rmeddis@38:                 % Separate MOC effect for each BF
rmeddis@38:                 HSRbegins=nBFs*(nANfiberTypes-1)+1;
rmeddis@38:                 rates=ICspikes(HSRbegins:end,:)/ANdt;
rmeddis@38:                 for idx=1:nBFs
rmeddis@38:                     [smoothedRates, MOCboundary{idx}] = ...
rmeddis@38:                         filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ...
rmeddis@38:                         MOCboundary{idx});
rmeddis@38:                     % spont 'rates' is zero for IC
rmeddis@38:                     MOCattSegment(idx,:)=smoothedRates;
rmeddis@38:                     % expand timescale back to model dt from ANdt
rmeddis@38:                     x= repmat(MOCattSegment(idx,:), ANspeedUpFactor,1);
rmeddis@38:                     x= reshape(x,1,segmentLength);
rmeddis@38:                     MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ...
rmeddis@38:                         (1- MOCrateToAttenuationFactor*  x);
rmeddis@38:                 end
rmeddis@38:                 MOCattenuation(MOCattenuation<0)=0.04;
rmeddis@38:                 % segment debugging
rmeddis@38:                 % plotInstructions.figureNo=98;
rmeddis@38:                 % plotInstructions.displaydt=ANdt;
rmeddis@38:                 %  plotInstructions.numPlots=1;
rmeddis@38:                 %  plotInstructions.subPlotNo=1;
rmeddis@38:                 % UTIL_plotMatrix(ICspikes, plotInstructions);
rmeddis@38: 
rmeddis@38:     end     % AN_spikesOrProbability
rmeddis@38:     segmentStartPTR=segmentStartPTR+segmentLength;
rmeddis@38:     reducedSegmentPTR=reducedSegmentPTR+reducedSegmentLength;
rmeddis@38: 
rmeddis@38: 
rmeddis@38: end  % segment
rmeddis@38: 
rmeddis@38: %% apply refractory correction to spike probabilities
rmeddis@38: 
rmeddis@38: % the probability of a spike's having occurred in the preceding 
rmeddis@38: %  refractory window 
rmeddis@38: %    pFired= 1 - II(1-p(t)), t= tnow-refractory period: tnow-1
rmeddis@38: % we need a running account of cumProb=II(1-p(t))
rmeddis@38: %   cumProb(t)= cumProb(t-1)*(1-p(t-1))/(1-p(t-refracPeriod))
rmeddis@38: %   cumProb(0)=0
rmeddis@38: %   pFired(t)= 1-cumProb(t)
rmeddis@38: % whence
rmeddis@38: %   p(t)= ANprobOutput(t) * pFired(t)
rmeddis@38: % where ANprobOutput is the uncorrected probability
rmeddis@38: 
rmeddis@38: 
rmeddis@38: % switch AN_spikesOrProbability
rmeddis@38: %     case 'probability'
rmeddis@38: %         ANprobOutput=ANprobRateOutput*dt;
rmeddis@38: %         [r nEpochs]=size(ANprobOutput);
rmeddis@38: %         % find probability of no spikes in refractory period
rmeddis@38: %         pNoSpikesInRefrac=ones(size(ANprobOutput));
rmeddis@38: %         pSpike=zeros(size(ANprobOutput));
rmeddis@38: %         for epochNo=lengthAbsRefractoryP+2:nEpochs
rmeddis@38: %             pNoSpikesInRefrac(:,epochNo)=...
rmeddis@38: %                 pNoSpikesInRefrac(:,epochNo-2)...
rmeddis@38: %                 .*(1-pSpike(:,epochNo-1))...
rmeddis@38: %                 ./(1-pSpike(:,epochNo-lengthAbsRefractoryP-1));
rmeddis@38: %             pSpike(:,epochNo)= ANprobOutput(:,epochNo)...
rmeddis@38: %                 .*pNoSpikesInRefrac(:,epochNo);
rmeddis@38: %         end
rmeddis@38: %         ANprobRateOutput=pSpike/dt;
rmeddis@38: % end
rmeddis@38: 
rmeddis@38: path(restorePath)