MAP

function [iih,IPIhisttime,IPIhistweight]=IPIHextract(IFRAN_pattern, sfreq)

% 

% tracks the formants according to an analysis proposed in Secker-Walker

% JASA 1990, section V.A

% Tim J?rgens, February 2011, code from Guy Brown included

% 

% input:  IFRAN_pattern: pattern of the auditory model (dependend on the number of modules used)

%                        first dimension: frequency channel, 

%                        second dimension: time (samples)

%         sfreq:         sampling frequency

% output: iih:           interpeak-interval histogram, matrix very similar

%                        the plot 5 in the Secker-Walker paper

%

%

%

time_axis = 0:1/sfreq:(size(IFRAN_pattern,2)-1)/sfreq;

%find how many samples of AN_pattern are 10ms and 3ms

%one_sample_is_a_time_of = time_axis(2);

[tmp, start_time_index] = min(abs(0-time_axis));

[tmp, stop20_time_index] = min(abs(0.020-time_axis));

number_of_samples20ms = stop20_time_index - start_time_index;

[tmp, stop10_time_index] = min(abs(0.010-time_axis));

number_of_samples10ms = stop10_time_index - start_time_index;

every_10ms = 1:number_of_samples10ms:size(IFRAN_pattern,2)-number_of_samples20ms;

hamm_window = hamming(11);

halfHamming = (length(hamm_window)-1)/2;

% window normalization

norm = conv(ones(1,number_of_samples20ms),hamm_window);

norm = norm(5+1:end-5)';

win_size = number_of_samples20ms;

half_win_size = floor(win_size/2);

hop_size = number_of_samples10ms;

%parameters of the autocorrelation

params.acfTau=0.1;

params.lags=[0:1/sfreq:0.02-1/sfreq];

sampledacf = runningACF(IFRAN_pattern,sfreq,params);

sampledacf(sampledacf<0)=0;

sampledacf = sqrt(sampledacf);

%pre-allocation due to speed

%Acorr = zeros(size(IFRAN_pattern,1),size(every_3ms,2),number_of_samples10ms*2+1);

%RAcorr = zeros(size(IFRAN_pattern,1),size(every_3ms,2),number_of_samples10ms*2+1);

%SRAcorr = zeros(size(IFRAN_pattern,1),size(every_3ms,2),number_of_samples10ms*2+1-10);

IPIhisttime = zeros(size(IFRAN_pattern,1),size(every_10ms,2),3);

IPIhistweight = zeros(size(IFRAN_pattern,1),size(every_10ms,2),3);  %maximum 3 peaks from the SRA

iih = zeros(half_win_size,size(sampledacf,1));

for iCounter = 1:size(sampledacf,2) %each channel

    fprintf('Channel No. %i\n',iCounter);

    %time_counter = 1;

    %for jCounter = every_3ms %every 3ms time segment

    for frame=1:size(sampledacf,1)

        %%debug

        %if iCounter == 130

        %    disp('here');

        %end

        sra = conv(squeeze(sampledacf(frame,iCounter,:)),hamm_window);

        sra = sra(halfHamming+1:end-halfHamming)./norm;

        df = [0 ; diff(sra)];

        idx = find((df(1:end-1)>=0)&(df(2:end)<0));

        % interpolate

        a=df(idx);

        b=df(idx+1);

        idx = (idx-1+a./(a-b));

        % get rid of a zero peak, if it exists

        idx = idx(idx>1);

        %include the zeroth peak

        idx = [1 idx']';

        % peak values corresponding to these intervals

        amp = interp1(1:length(sra),sra,idx,'linear');

        % if required, remove peaks that lie below the mean sra

        % note that we disregard the value at zero delay

        %if (params.removePeaksBelowMean)

        valid = find(amp>1.2*mean(sra(1:floor(length(sra)/2))));

        %valid = find(amp>mean(sra));

        %just take the mean of the first half of the sra as a comparison

        idx = idx(valid);

        amp = amp(valid);

        %end

        % only use the first four peaks (three intervals)

        idx = idx(1:min(4,length(idx)));

        % find the intervals

        interval = diff(idx);

        % now histogram the intervals

        if (~isempty(interval))

            for k=1:length(interval)

                if interval(k)<=half_win_size

                    iih(round(interval(k)),frame) = iih(round(interval(k)),frame)+amp(k);

                    IPIhisttime(iCounter,frame,k) = interval(k)/sfreq;

                    IPIhistweight(iCounter,frame,k) = amp(k);

                end

            end

        end

    end

    %% end Guy's code

    %         %take the autocorrelation (ACF) of a 10ms-segment of each channel

    %         Acorr(iCounter,time_counter,:) = xcorr(IFRAN_pattern(iCounter,jCounter:number_of_samples10ms+jCounter),'biased'); %biased scales the ACF by the reciprocal of the length of the segment

    %         %root calculation

    %         RAcorr(iCounter,time_counter,:) = sqrt(abs(Acorr(iCounter,time_counter,:)));

    %

    %         %smoothing using the 11-point hamming window

    %         for kCounter = 6:size(RAcorr(iCounter,time_counter,:),3)-5 %start with 6 and end with 5 samples

    %             %less the length of time_axis not to get in conflict with the length of

    %             %the hamm_window

    %             SRAcorr(iCounter,time_counter,kCounter-5) = ...

    %                 squeeze(RAcorr(iCounter,time_counter,(kCounter-5):(kCounter+5)))'*hamm_window./sum(hamm_window);

    %         end

    %

    %         %mean value of actual SRA

    %         SRA_mean = mean(SRAcorr(iCounter,time_counter,:));

    %

    %         %find signed zero-crossings of the first derivative (=difference)

    %         z_crossings_indices = find(diff(sign(diff(squeeze(SRAcorr(iCounter,time_counter,:))))) < 0)+1; %+1 is necessary, because diff shortens vector by 1

    %         middle_index = ceil(size(SRAcorr(iCounter,time_counter,:),3)/2);

    %

    %         validCounter = 1;

    %         valid_z_crossings_indices = [];

    %         %find valid zero-crossings (peak higher than meanvalue and within first 5 ms of SRA)

    %         for lCounter = 1:length(z_crossings_indices)

    %             if (SRAcorr(iCounter,time_counter,z_crossings_indices(lCounter)) > SRA_mean) && ...

    %                     (abs(z_crossings_indices(lCounter)-middle_index) < round(number_of_samples10ms/2));

    %                 valid_z_crossings_indices(validCounter) = z_crossings_indices(lCounter);

    %                 validCounter = validCounter+1;

    %             end

    %         end

    %

    %         %find main peak in the ACF

    %         [tmp,index_of_z_crossings_main_index] = min(abs(middle_index-valid_z_crossings_indices));

    %         if ~tmp == 0

    %             disp('middle peak not appropriately found');

    %         end

    %

    %         %%% for debugging

    % %         if iCounter == 130

    % %             disp('here');

    % %             figure, plot(squeeze(SRAcorr(iCounter,time_counter,:)));

    % %             hold on, plot([1 length(squeeze(SRAcorr(iCounter,time_counter,:)))],[SRA_mean SRA_mean],'r-');

    % %         end

    %         %%%

    %

    %         %generate IPI-histogram: take the first 3 intervals of SRAcorr

    %         %(positive delay) in the first 5 ms

    %         histcounter = 1;

    %         for lCounter = index_of_z_crossings_main_index+1:min([length(valid_z_crossings_indices(index_of_z_crossings_main_index+1:end)) 3])+index_of_z_crossings_main_index

    %             sampledifference = abs(valid_z_crossings_indices(lCounter)-valid_z_crossings_indices(lCounter-1));

    %             %the difference between two adjacent peaks in the SRA is taken

    %             %as IPI estimate

    %             IPIhisttime(iCounter,time_counter,histcounter) = sampledifference*one_sample_is_a_time_of;

    %             %the amplitude of the SRA at the start of the SRA interval is

    %             %taken as the IPIweight

    %             IPIhistweight(iCounter,time_counter,histcounter) = SRAcorr(iCounter,time_counter,valid_z_crossings_indices(lCounter-1));

    %             histcounter = histcounter + 1;

    %         end

    %time_counter = time_counter+1;

end

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)

function Pavel_MAP1_14

% test_MAP1_14 is a general purpose test routine that can be adjusted to

% test a number of different applications of MAP1_14

%

% A range of options are supplied in the early part of the program

%

% One use of the function is to create demonstrations; filenames <demoxx>

%  to illustrate particular features

%

% #1

% Identify the file (in 'MAPparamsName') containing the model parameters

%

% #2

% Identify the kind of model required (in 'AN_spikesOrProbability').

%  A full brainstem model (spikes) can be computed or a shorter model

%  (probability) that computes only so far as the auditory nerve

%

% #3

% Choose between a tone signal or file input (in 'signalType')

%

% #4

% Set the signal rms level (in leveldBSPL)

%

% #5

% Identify the channels in terms of their best frequencies in the vector

%  BFlist.

%

% Last minute changes to the parameters fetched earlier can be made using

%  the cell array of strings 'paramChanges'.

%  Each string must have the same format as the corresponding line in the

%  file identified in 'MAPparamsName'

%

% When the demonstration is satisfactory, freeze it by renaming it <demoxx>

restorePath=path;

addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'], ...

    ['..' filesep 'utilities'])

%%  #1 parameter file name