MAP

function  MAP1_14(inputSignal, sampleRate, BFlist, MAPparamsName, ...

    AN_spikesOrProbability, paramChanges)

% To test this function use test_MAP1_14 in this folder

%

% All arguments are mandatory.

%

% BFlist is a list of BFs but can be '-1' to allow MAPparams to choose

%

%  MAPparamsName='Normal'; % source of model parameters

%  AN_spikesOrProbability='spikes'; % or 'probability'

%  paramChanges is a cell array of strings that can be used to make last

%   minute parameter changes, e.g., to simulate OHC loss

%   paramChanges{1}= 'DRNLParams.a=0;';

% The model parameters are established in the MAPparams<***> file

%  and stored as global

restorePath=path;

addpath (['..' filesep 'parameterStore'])

global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams

global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams

% All of the results of this function are stored as global

global dt ANdt  savedBFlist saveAN_spikesOrProbability saveMAPparamsName...

    savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ...

    DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...

    IHCoutput ANprobRateOutput ANoutput savePavailable tauCas  ...

    CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation

% Normally only ICoutput(logical spike matrix) or ANprobRateOutput will be

% needed by the user; so the following will suffice

%   global ANdt ICoutput ANprobRateOutput

% Note that sampleRate has not changed from the original function call and

%  ANprobRateOutput is sampled at this rate

% However ANoutput, CNoutput and IC output are stored as logical

%  'spike' matrices using a lower sample rate (see ANdt).

% When AN_spikesOrProbability is set to probability,

%  no spike matrices are computed.

% When AN_spikesOrProbability is set to 'spikes',

%  no probability output is computed

% Efferent control variables are ARattenuation and MOCattenuation

%  These are scalars between 1 (no attenuation) and 0.

%  They are represented with dt=1/sampleRate (not ANdt)

%  They are computed using either AN probability rate output

%   or IC (spikes) output as approrpriate.

% AR is computed using across channel activity

% MOC is computed on a within-channel basis.

% save as global for later plotting if required

savedBFlist=BFlist;

saveAN_spikesOrProbability=AN_spikesOrProbability;

saveMAPparamsName=MAPparamsName;

% Read parameters from MAPparams<***> file in 'parameterStore' folder

cmd=['method=MAPparams' MAPparamsName ...

    '(BFlist, sampleRate, 0);'];

eval(cmd);

% Beware, 'BFlist=-1' is a legitimate argument for MAPparams<>

%  if the calling program allows MAPparams to specify the list

BFlist=DRNLParams.nonlinCFs;

% now accept last mintue parameter changes required by the calling program

if nargin>5 && ~isempty(paramChanges)

    nChanges=length(paramChanges);

    for idx=1:nChanges

        eval(paramChanges{idx})

    end

end

dt=1/sampleRate;

duration=length(inputSignal)/sampleRate;

% segmentDuration is specified in parameter file (must be >efferent delay)

segmentDuration=method.segmentDuration;

segmentLength=round(segmentDuration/ dt);

segmentTime=dt*(1:segmentLength); % used in debugging plots

% all spiking activity is computed using longer epochs

ANspeedUpFactor=5;  % 5 times longer

% 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

% Nyquist=(1/ANdt)/2;

% [ARfilt_b,ARfilt_a] = butter(1, (1/(2*pi*OMEParams.ARtau))/Nyquist, 'low');

a1=dt/OMEParams.ARtau-1; a0=1;

b0=1+ a1;

ARfilt_b=b0; ARfilt_a=[a0 a1];

ARattenuation=ones(1,signalLength);

ARrateThreshold=OMEParams.ARrateThreshold; % may not be used

ARrateToAttenuationFactor=OMEParams.rateToAttenuationFactor;

ARrateToAttenuationFactorProb=OMEParams.rateToAttenuationFactorProb;

ARboundary=[];

ARboundaryProb=0;

% save complete OME record (stapes displacement)

OMEoutput=zeros(1,signalLength);

TMoutput=zeros(1,signalLength);

%% BM ---

% BM is represented as a list of locations identified by BF

DRNL_BFs=BFlist;

nBFs= length(DRNL_BFs);

% DRNLchannelParameters=DRNLParams.channelParameters;

DRNLresponse= zeros(nBFs, segmentLength);

MOCrateToAttenuationFactor=DRNLParams.rateToAttenuationFactor;

rateToAttenuationFactorProb=DRNLParams.rateToAttenuationFactorProb;

MOCrateThreshold=DRNLParams.MOCrateThreshold;

% smoothing filter for MOC

% Nyquist=(1/ANdt)/2;

% [MOCfilt_b,MOCfilt_a] = ...

%     butter(1, (1/(2*pi*DRNLParams.MOCtau))/Nyquist, 'low');

% figure(10), freqz(stapesDisp_b, stapesDisp_a)

a1=dt/DRNLParams.MOCtau-1; a0=1;

b0=1+ a1;

MOCfilt_b=b0; MOCfilt_a=[a0 a1];

% figure(9), freqz(stapesDisp_b, stapesDisp_a)

MOCboundary=cell(nBFs,1);

MOCprobBoundary=cell(nBFs,1);

MOCattSegment=zeros(nBFs,reducedSegmentLength);

MOCattenuation=ones(nBFs,signalLength);

if DRNLParams.a>0

    DRNLcompressionThreshold=10^((1/(1-DRNLParams.c))* ...

    log10(DRNLParams.b/DRNLParams.a));

else

    DRNLcompressionThreshold=inf;

end

DRNLlinearOrder= DRNLParams.linOrder;

DRNLnonlinearOrder= DRNLParams.nonlinOrder;

DRNLa=DRNLParams.a;

DRNLb=DRNLParams.b;

DRNLc=DRNLParams.c;

linGAIN=DRNLParams.g;

%

% gammatone filter coefficients for linear pathway

bw=DRNLParams.linBWs';

phi = 2 * pi * bw * dt;

cf=DRNLParams.linCFs';

theta = 2 * pi * cf * dt;

cos_theta = cos(theta);

sin_theta = sin(theta);

alpha = -exp(-phi).* cos_theta;

b0 = ones(nBFs,1);

b1 = 2 * alpha;

b2 = exp(-2 * phi);

z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i;

z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i;

z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i;

tf = (z2 + z3) ./ z1;

a0 = abs(tf);

a1 = alpha .* a0;

GTlin_a = [b0, b1, b2];

GTlin_b = [a0, a1];

GTlinOrder=DRNLlinearOrder;

GTlinBdry=cell(nBFs,GTlinOrder);

% nonlinear gammatone filter coefficients 

bw=DRNLParams.nlBWs';

phi = 2 * pi * bw * dt;

cf=DRNLParams.nonlinCFs';

theta = 2 * pi * cf * dt;

cos_theta = cos(theta);

sin_theta = sin(theta);

alpha = -exp(-phi).* cos_theta;

b0 = ones(nBFs,1);

b1 = 2 * alpha;

b2 = exp(-2 * phi);

z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i;

z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i;

z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i;

tf = (z2 + z3) ./ z1;

a0 = abs(tf);

a1 = alpha .* a0;

GTnonlin_a = [b0, b1, b2];

GTnonlin_b = [a0, a1];

GTnonlinOrder=DRNLnonlinearOrder;

GTnonlinBdry1=cell(nBFs, GTnonlinOrder);

GTnonlinBdry2=cell(nBFs, GTnonlinOrder);

% complete BM record (BM displacement)

DRNLoutput=zeros(nBFs, signalLength);

%% IHC ---

% IHC cilia activity and receptor potential

% viscous coupling between BM and stereocilia displacement

% Nyquist=sampleRate/2;

% IHCcutoff=1/(2*pi*IHC_cilia_RPParams.tc);

% [IHCciliaFilter_b,IHCciliaFilter_a]=...

%     butter(1, IHCcutoff/Nyquist, 'high');

a1=dt/IHC_cilia_RPParams.tc-1; a0=1;

b0=1+ a1;

% high pass (i.e. low pass reversed)

IHCciliaFilter_b=[a0 a1]; IHCciliaFilter_a=b0;

% figure(9), freqz(IHCciliaFilter_b, IHCciliaFilter_a)

IHCciliaBndry=cell(nBFs,1);

% IHC apical conductance (Boltzman function)

IHC_C= IHC_cilia_RPParams.C;

IHCu0= IHC_cilia_RPParams.u0;

IHCu1= IHC_cilia_RPParams.u1;

IHCs0= IHC_cilia_RPParams.s0;

IHCs1= IHC_cilia_RPParams.s1;

IHCGmax= IHC_cilia_RPParams.Gmax;

IHCGu0= IHC_cilia_RPParams.Gu0; % (leakage)

IHCGa= IHCGmax./(1+exp(-(0-IHCu0)/IHCs0).*(1+exp(-(0-IHCu1)/IHCs1)));

IHCrestingCiliaCond=IHCGa+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= -0.06;

IHC_Vnow= IHCrestingV*ones(nBFs,1); % initial voltage

IHC_RP= zeros(nBFs,segmentLength);

% complete record of IHC receptor potential (V)

IHCciliaDisplacement= zeros(nBFs,segmentLength);

IHCoutput= zeros(nBFs,signalLength);

IHC_cilia_output= zeros(nBFs,signalLength);

%% pre-synapse ---

% Each BF is replicated using a different fiber type to make a 'channel'

% The number of channels is nBFs x nANfiberTypes

% Fiber types are specified in terms of tauCa

nANfiberTypes= length(IHCpreSynapseParams.tauCa);

tauCas= IHCpreSynapseParams.tauCa;

nChannels= nANfiberTypes*nBFs;

synapticCa= zeros(nChannels,segmentLength);

% Calcium control (more calcium, greater release rate)

ECa=IHCpreSynapseParams.ECa;

gamma=IHCpreSynapseParams.gamma;

beta=IHCpreSynapseParams.beta;

tauM=IHCpreSynapseParams.tauM;

mICa=zeros(nChannels,segmentLength);

GmaxCa=IHCpreSynapseParams.GmaxCa;

synapse_z= IHCpreSynapseParams.z;

synapse_power=IHCpreSynapseParams.power;

% tauCa vector is established across channels to allow vectorization

%  (one tauCa per channel). Do not confuse with tauCas (one pre fiber type)

tauCa=repmat(tauCas, nBFs,1);

tauCa=reshape(tauCa, nChannels, 1);

% presynapse startup values (vectors, length:nChannels)

% proportion (0 - 1) of Ca channels open at IHCrestingV

mICaCurrent=((1+beta^-1 * exp(-gamma*IHCrestingV))^-1)...

    *ones(nBFs*nANfiberTypes,1);

% corresponding startup currents

ICaCurrent= (GmaxCa*mICaCurrent.^3) * (IHCrestingV-ECa);

CaCurrent= ICaCurrent.*tauCa;

% vesicle release rate at startup (one per channel)

% kt0 is used only at initialisation

kt0= -synapse_z * CaCurrent.^synapse_power;

%% AN ---

% each row of the AN matrices represents one AN fiber

% The results computed either for probabiities *or* for spikes (not both)

% Spikes are necessary if CN and IC are to be computed

nFibersPerChannel= AN_IHCsynapseParams.numFibers;

nANfibers= nChannels*nFibersPerChannel;

y=AN_IHCsynapseParams.y;

l=AN_IHCsynapseParams.l;

x=AN_IHCsynapseParams.x;

r=AN_IHCsynapseParams.r;

M=round(AN_IHCsynapseParams.M);

% probability            (NB initial 'P' on everything)

PAN_ydt = repmat(AN_IHCsynapseParams.y*dt, nChannels,1);

PAN_ldt = repmat(AN_IHCsynapseParams.l*dt, nChannels,1);

PAN_xdt = repmat(AN_IHCsynapseParams.x*dt, nChannels,1);

PAN_rdt = repmat(AN_IHCsynapseParams.r*dt, nChannels,1);

PAN_rdt_plus_ldt = PAN_rdt + PAN_ldt;

PAN_M=round(AN_IHCsynapseParams.M);

% compute starting values

Pcleft    = kt0* y* M ./ (y*(l+r)+ kt0* l);

Pavailable    = Pcleft*(l+r)./kt0;

Preprocess    = Pcleft*r/x; % canbe fractional

ANprobability=zeros(nChannels,segmentLength);

ANprobRateOutput=zeros(nChannels,signalLength);

% special variables for monitoring synaptic cleft (specialists only)

savePavailableSeg=zeros(nChannels,segmentLength);

savePavailable=zeros(nChannels,signalLength);

% spikes     % !  !  !    ! !        !   !  !

AN_refractory_period= AN_IHCsynapseParams.refractory_period;

lengthAbsRefractory= round(AN_refractory_period/ANdt);

AN_ydt= repmat(AN_IHCsynapseParams.y*ANdt, nANfibers,1);

AN_ldt= repmat(AN_IHCsynapseParams.l*ANdt, nANfibers,1);

AN_xdt= repmat(AN_IHCsynapseParams.x*ANdt, nANfibers,1);

AN_rdt= repmat(AN_IHCsynapseParams.r*ANdt, nANfibers,1);

AN_rdt_plus_ldt= AN_rdt + AN_ldt;

AN_M= round(AN_IHCsynapseParams.M);

% kt0  is initial release rate

% Establish as a vector (length=channel x number of fibers)

kt0= repmat(kt0', nFibersPerChannel, 1);

kt0=reshape(kt0, nANfibers,1);

% starting values for reservoirs

AN_cleft    = kt0* y* M ./ (y*(l+r)+ kt0* l);

AN_available    = round(AN_cleft*(l+r)./kt0); %must be integer

AN_reprocess    = AN_cleft*r/x;

% output is in a logical array spikes = 1/ 0.

ANspikes= false(nANfibers,reducedSegmentLength);

ANoutput= false(nANfibers,reducedSignalLength);

%% CN (first brain stem nucleus - could be any subdivision of CN)

% Input to a CN neuorn is a random selection of AN fibers within a channel

%  The number of AN fibers used is ANfibersFanInToCN

ANfibersFanInToCN=MacGregorMultiParams.fibersPerNeuron;

nCNneuronsPerChannel=MacGregorMultiParams.nNeuronsPerBF;

% CNtauGk (Potassium time constant) determines the rate of firing of

%  the unit when driven hard by a DC input (not normally >350 sp/s)

CNtauGk=MacGregorMultiParams.tauGk;

ANavailableFibersPerChan=AN_IHCsynapseParams.numFibers;

nCNneurons=nCNneuronsPerChannel*nChannels;

% nCNneuronsPerFiberType= nCNneurons/nANfiberTypes;

CNmembranePotential=zeros(nCNneurons,reducedSegmentLength);

% establish which ANfibers (by name) feed into which CN nuerons

CNinputfiberLists=zeros(nChannels*nCNneuronsPerChannel, ANfibersFanInToCN);

unitNo=1;

for ch=1:nChannels

    % Each channel contains a number of units =length(listOfFanInValues)

    for idx=1:nCNneuronsPerChannel

        fibersUsed=(ch-1)*ANavailableFibersPerChan + ...

            ceil(rand(1,ANfibersFanInToCN)* ANavailableFibersPerChan);

        CNinputfiberLists(unitNo,:)=fibersUsed;

        unitNo=unitNo+1;

    end

end

% input to CN units

AN_PSTH=zeros(nCNneurons,reducedSegmentLength);

% Generate CNalphaFunction function

%  by which spikes are converted to post-synaptic currents

CNdendriteLPfreq= MacGregorMultiParams.dendriteLPfreq;

CNcurrentPerSpike=MacGregorMultiParams.currentPerSpike;

CNspikeToCurrentTau=1/(2*pi*CNdendriteLPfreq);

t=ANdt:ANdt:5*CNspikeToCurrentTau;

CNalphaFunction=...

    (CNcurrentPerSpike/CNspikeToCurrentTau)*t.*exp(-t/CNspikeToCurrentTau);

% figure(98), plot(t,CNalphaFunction)

% working memory for implementing convolution

CNcurrentTemp=...

    zeros(nCNneurons,reducedSegmentLength+length(CNalphaFunction)-1);

% trailing alphas are parts of humps carried forward to the next segment

CNtrailingAlphas=zeros(nCNneurons,length(CNalphaFunction));

CN_tauM=MacGregorMultiParams.tauM;

CN_tauTh=MacGregorMultiParams.tauTh;

CN_cap=MacGregorMultiParams.Cap;

CN_c=MacGregorMultiParams.c;

CN_b=MacGregorMultiParams.dGkSpike;

CN_Ek=MacGregorMultiParams.Ek;

CN_Eb= MacGregorMultiParams.Eb;

CN_Er=MacGregorMultiParams.Er;

CN_Th0= MacGregorMultiParams.Th0;

CN_E= zeros(nCNneurons,1);

CN_Gk= zeros(nCNneurons,1);

CN_Th= MacGregorMultiParams.Th0*ones(nCNneurons,1);

CN_Eb=CN_Eb.*ones(nCNneurons,1);

CN_Er=CN_Er.*ones(nCNneurons,1);

CNtimeSinceLastSpike=zeros(nCNneurons,1);

% tauGk is the main distinction between neurons

%  in fact they are all the same in the standard model

tauGk=repmat(CNtauGk,nChannels*nCNneuronsPerChannel,1);

CN_PSTH=zeros(nChannels,reducedSegmentLength);

CNoutput=false(nCNneurons,reducedSignalLength);

%% MacGregor (IC - second nucleus) --------

nICcells=nChannels;  % one cell per channel

ICspikeWidth=0.00015;   % this may need revisiting

epochsPerSpike=round(ICspikeWidth/ANdt);

if epochsPerSpike<1

    error(['MacGregorMulti: sample rate too low to support ' ...

        num2str(ICspikeWidth*1e6) '  microsec spikes']);

end

% short names

IC_tauM=MacGregorParams.tauM;

IC_tauGk=MacGregorParams.tauGk;

IC_tauTh=MacGregorParams.tauTh;

IC_cap=MacGregorParams.Cap;

IC_c=MacGregorParams.c;

IC_b=MacGregorParams.dGkSpike;

IC_Th0=MacGregorParams.Th0;

IC_Ek=MacGregorParams.Ek;

IC_Eb= MacGregorParams.Eb;

IC_Er=MacGregorParams.Er;

IC_E=zeros(nICcells,1);

IC_Gk=zeros(nICcells,1);

IC_Th=IC_Th0*ones(nICcells,1);

% Dendritic filtering, all spikes are replaced by CNalphaFunction functions

ICdendriteLPfreq= MacGregorParams.dendriteLPfreq;

ICcurrentPerSpike=MacGregorParams.currentPerSpike;

ICspikeToCurrentTau=1/(2*pi*ICdendriteLPfreq);

t=ANdt:ANdt:3*ICspikeToCurrentTau;

IC_CNalphaFunction= (ICcurrentPerSpike / ...

    ICspikeToCurrentTau)*t.*exp(-t / ICspikeToCurrentTau);

% figure(98), plot(t,IC_CNalphaFunction)

% working space for implementing alpha function

ICcurrentTemp=...

    zeros(nICcells,reducedSegmentLength+length(IC_CNalphaFunction)-1);

ICtrailingAlphas=zeros(nICcells, length(IC_CNalphaFunction));

ICfiberTypeRates=zeros(nANfiberTypes,reducedSignalLength);

ICoutput=false(nChannels,reducedSignalLength);

ICmembranePotential=zeros(nICcells,reducedSegmentLength);

ICmembraneOutput=zeros(nICcells,signalLength);

%% Main program %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%

%  Compute the entire model for each segment

segmentStartPTR=1;

reducedSegmentPTR=1; % when sampling rate is reduced

while segmentStartPTR<signalLength

    segmentEndPTR=segmentStartPTR+segmentLength-1;

    % shorter segments after speed up.

    shorterSegmentEndPTR=reducedSegmentPTR+reducedSegmentLength-1;

    iputPressureSegment=inputSignal...

        (:,segmentStartPTR:segmentStartPTR+segmentLength-1);

    % segment debugging plots

    % figure(98)

    % plot(segmentTime,iputPressureSegment), title('signalSegment')

    % OME ----------------------

    % OME Stage 1: external resonances. Add to inputSignal pressure wave

    y=iputPressureSegment;

    for n=1:nOMEExtFilters

        % any number of resonances can be used

        [x  OMEExtFilterBndry{n}] = ...

            filter(ExtFilter_b{n},ExtFilter_a{n},...

            iputPressureSegment, OMEExtFilterBndry{n});

        x= x* OMEgainScalars(n);

        % This is a parallel resonance so add it

        y=y+x;

    end

    iputPressureSegment=y;

    OMEextEarPressure(segmentStartPTR:segmentEndPTR)= iputPressureSegment;

    % OME stage 2: convert input pressure (velocity) to

    %  tympanic membrane(TM) displacement using low pass filter

    [TMdisplacementSegment  OME_TMdisplacementBndry] = ...

        filter(TMdisp_b,TMdisp_a,iputPressureSegment, ...

        OME_TMdisplacementBndry);

    % and save it

    TMoutput(segmentStartPTR:segmentEndPTR)= TMdisplacementSegment;

    % OME stage 3: middle ear high pass effect to simulate stapes inertia

    [stapesDisplacement  OMEhighPassBndry] = ...

        filter(stapesDisp_b,stapesDisp_a,TMdisplacementSegment, ...

        OMEhighPassBndry);

    % OME stage 4:  apply stapes scalar

    stapesDisplacement=stapesDisplacement*stapesScalar;

    % OME stage 5:    acoustic reflex stapes attenuation

    %  Attenuate the TM response using feedback from LSR fiber activity

    if segmentStartPTR>efferentDelayPts

        stapesDisplacement= stapesDisplacement.*...

            ARattenuation(segmentStartPTR-efferentDelayPts:...

            segmentEndPTR-efferentDelayPts);

    end

    % segment debugging plots

    % figure(98)

    % plot(segmentTime, stapesDisplacement), title ('stapesDisplacement')

    % and save

    OMEoutput(segmentStartPTR:segmentEndPTR)= stapesDisplacement;

    %% BM ------------------------------

    % Each location is computed separately

    for BFno=1:nBFs

        %            *linear* path

        linOutput = stapesDisplacement * linGAIN;  % linear gain

        for order = 1 : GTlinOrder

            [linOutput GTlinBdry{BFno,order}] = ...

                filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, GTlinBdry{BFno,order});

        end

        %           *nonLinear* path

        % efferent attenuation (0 <> 1)

        if segmentStartPTR>efferentDelayPts

            MOC=MOCattenuation(BFno, segmentStartPTR-efferentDelayPts:...

                segmentEndPTR-efferentDelayPts);

        else    % no MOC available yet

            MOC=ones(1, segmentLength);

        end

        %       first gammatone filter

        for order = 1 : GTnonlinOrder

            [nonlinOutput GTnonlinBdry1{BFno,order}] = ...

                filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ...

                stapesDisplacement, GTnonlinBdry1{BFno,order});

        end

        %       broken stick instantaneous compression

        % nonlinear gain is weakend by MOC before applied to BM response

        y= nonlinOutput.*(MOC* DRNLa);  % linear section.

        % compress those parts of the signal above the compression

        % threshold

        abs_x = abs(nonlinOutput);

        idx=find(abs_x>DRNLcompressionThreshold);

        if ~isempty(idx)>0

            y(idx)=sign(nonlinOutput(idx)).*...

                (DRNLb*abs_x(idx).^DRNLc);

        end

        nonlinOutput=y;

        %       second filter removes distortion products

        for order = 1 : GTnonlinOrder

            [ nonlinOutput GTnonlinBdry2{BFno,order}] = ...

                filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), nonlinOutput, GTnonlinBdry2{BFno,order});

        end

        %  combine the two paths to give the DRNL displacement

        DRNLresponse(BFno,:)=linOutput+nonlinOutput;

    end % BF

    % segment debugging plots

    % figure(98)

    %     if size(DRNLresponse,1)>3

    %         imagesc(DRNLresponse)  % matrix display

    %         title('DRNLresponse'); % single or double channel response

    %     else

    %         plot(segmentTime, DRNLresponse)

    %     end

    % and save it

    DRNLoutput(:, segmentStartPTR:segmentEndPTR)= DRNLresponse;

    %% IHC ------------------------------------

    %  BM displacement to IHCciliaDisplacement is  a high-pass filter

    %   because of viscous coupling

    for idx=1:nBFs

        [IHCciliaDisplacement(idx,:)  IHCciliaBndry{idx}] = ...

            filter(IHCciliaFilter_b,IHCciliaFilter_a, ...

            DRNLresponse(idx,:), IHCciliaBndry{idx});

    end

    % apply scalar

    IHCciliaDisplacement=IHCciliaDisplacement* IHC_C;

    % and save it

    IHC_cilia_output(:,segmentStartPTR:segmentStartPTR+segmentLength-1)=...

        IHCciliaDisplacement;

    % compute apical conductance

    G=1./(1+exp(-(IHCciliaDisplacement-IHCu0)/IHCs0).*...

        (1+exp(-(IHCciliaDisplacement-IHCu1)/IHCs1)));

    Gu=IHCGmax*G;

    % add resting conductance to give apical conductance

    Gu= Gu+IHCGu0;

    % Compute receptor potential

    for idx=1:segmentLength

        IHC_Vnow=IHC_Vnow+ (-Gu(:, idx).*(IHC_Vnow-IHC_Et)-...

            IHC_Gk*(IHC_Vnow-IHC_Ekp))*  dt/IHC_Cab;

        IHC_RP(:,idx)=IHC_Vnow;

    end

    % segment debugging plots

    %     if size(IHC_RP,1)>3

    %         surf(IHC_RP), shading interp, title('IHC_RP')

    %     else

    %         plot(segmentTime, IHC_RP)

    %     end

    % and save it

    IHCoutput(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=IHC_RP;

    %% synapse -----------------------------

    % Compute the vesicle release rate for each fiber type at each BF

    % replicate IHC_RP for each fiber type

    Vsynapse=repmat(IHC_RP, nANfiberTypes,1);

    % look-up table of target fraction channels open for a given IHC_RP

    mICaINF=    1./( 1 + exp(-gamma  * Vsynapse)  /beta);

    % fraction of channel open - apply time constant

    for idx=1:segmentLength

        % mICaINF is the current 'target' value of mICa

        mICaCurrent=mICaCurrent+(mICaINF(:,idx)-mICaCurrent)*dt./tauM;

        mICa(:,idx)=mICaCurrent;

    end

    ICa=   (GmaxCa* mICa.^3) .* (Vsynapse- ECa);

    for idx=1:segmentLength

        CaCurrent=CaCurrent +  ICa(:,idx)*dt - CaCurrent*dt./tauCa;

        synapticCa(:,idx)=CaCurrent;

    end

    synapticCa=-synapticCa; % treat IHCpreSynapseParams as positive substance

    % NB vesicleReleaseRate is /s and is independent of dt

    vesicleReleaseRate = synapse_z * synapticCa.^synapse_power; % rate

    % segment debugging plots

    %     if size(vesicleReleaseRate,1)>3

    %         surf(vesicleReleaseRate), shading interp, title('vesicleReleaseRate')

    %     else

    %         plot(segmentTime, vesicleReleaseRate)

    %     end

    %% AN

    switch AN_spikesOrProbability

        case 'probability'

            % No refractory effect is applied

            for t = 1:segmentLength;

                M_Pq=PAN_M-Pavailable;

                M_Pq(M_Pq<0)=0;

                Preplenish = M_Pq .* PAN_ydt;

                Pejected = Pavailable.* vesicleReleaseRate(:,t)*dt;

                Preprocessed = M_Pq.*Preprocess.* PAN_xdt;

                ANprobability(:,t)= min(Pejected,1);

                reuptakeandlost= PAN_rdt_plus_ldt .* Pcleft;

                reuptake= PAN_rdt.* Pcleft;

                Pavailable= Pavailable+ Preplenish- Pejected+ Preprocessed;

                Pcleft= Pcleft + Pejected - reuptakeandlost;

                Preprocess= Preprocess + reuptake - Preprocessed;

                Pavailable(Pavailable<0)=0;

                savePavailableSeg(:,t)=Pavailable;    % synapse tracking

            end

            % and save it as *rate*

            ANrate=ANprobability/dt;

            ANprobRateOutput(:, segmentStartPTR:...

                segmentStartPTR+segmentLength-1)=  ANrate;

            % monitor synapse contents (only sometimes used)

            savePavailable(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=...

                savePavailableSeg;

            % Estimate efferent effects. ARattenuation (0 <> 1)

            %  acoustic reflex

            ARAttSeg=mean(ANrate(1:nBFs,:),1); %LSR channels are 1:nBF

            % smooth

            [ARAttSeg, ARboundaryProb] = ...

                filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundaryProb);

            ARAttSeg=ARAttSeg-ARrateThreshold;

            ARAttSeg(ARAttSeg<0)=0;   % prevent negative strengths

            ARattenuation(segmentStartPTR:segmentEndPTR)=...

                (1-ARrateToAttenuationFactorProb.* ARAttSeg);

            % MOC attenuation

            % within-channel HSR response only

            HSRbegins=nBFs*(nANfiberTypes-1)+1;

            rates=ANrate(HSRbegins:end,:);

            for idx=1:nBFs

                [smoothedRates, MOCprobBoundary{idx}] = ...

                    filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ...

                    MOCprobBoundary{idx});

                smoothedRates=smoothedRates-MOCrateThreshold;

                smoothedRates(smoothedRates<0)=0;

                MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ...

                    (1- smoothedRates* rateToAttenuationFactorProb);

            end

            MOCattenuation(MOCattenuation<0)=0.001;

        case 'spikes'

            ANtimeCount=0;

            % implement speed upt

            for t = ANspeedUpFactor:ANspeedUpFactor:segmentLength;

                ANtimeCount=ANtimeCount+1;

                % convert release rate to probabilities

                releaseProb=vesicleReleaseRate(:,t)*ANdt;

                % releaseProb is the release probability per channel

                %  but each channel has many synapses

                releaseProb=repmat(releaseProb',nFibersPerChannel,1);

                releaseProb=reshape(releaseProb, nFibersPerChannel*nChannels,1);

                % AN_available=round(AN_available); % vesicles must be integer, (?needed)

                M_q=AN_M- AN_available;     % number of missing vesicles

                M_q(M_q<0)= 0;              % cannot be less than 0

                % AN_N1 converts probability to discrete events

                %   it considers each event that might occur

                %   (how many vesicles might be released)

                %   and returns a count of how many were released

                % slow line

%                 probabilities= 1-(1-releaseProb).^AN_available;

                probabilities= 1-intpow((1-releaseProb), AN_available);

                ejected= probabilities> rand(length(AN_available),1);

                reuptakeandlost = AN_rdt_plus_ldt .* AN_cleft;

                reuptake = AN_rdt.* AN_cleft;

                % slow line

%                 probabilities= 1-(1-AN_reprocess.*AN_xdt).^M_q;

                probabilities= 1-intpow((1-AN_reprocess.*AN_xdt), M_q);

                reprocessed= probabilities>rand(length(M_q),1);

                % slow line

%                 probabilities= 1-(1-AN_ydt).^M_q;

                 probabilities= 1-intpow((1-AN_ydt), M_q);

                replenish= probabilities>rand(length(M_q),1);

                AN_available = AN_available + replenish - ejected ...

                    + reprocessed;

                AN_cleft = AN_cleft + ejected - reuptakeandlost;

                AN_reprocess = AN_reprocess + reuptake - reprocessed;

                % ANspikes is logical record of vesicle release events>0

                ANspikes(:, ANtimeCount)= ejected;

            end % t

            % zero any events that are preceded by release events ...

            %  within the refractory period

            % The refractory period consist of two periods

            %  1) the absolute period where no spikes occur

            %  2) a relative period where a spike may occur. This relative

            %     period is realised as a variable length interval

            %     where the length is chosen at random

            %     (esentially a linear ramp up)

            % Andreas has a fix for this

            for t = 1:ANtimeCount-2*lengthAbsRefractory;

                % identify all spikes across fiber array at time (t)

                % idx is a list of channels where spikes occurred

                % ?? try sparse matrices?

                idx=find(ANspikes(:,t));

                for j=idx  % consider each spike

                    % specify variable refractory period

                    %  between abs and 2*abs refractory period

                    nPointsRefractory=lengthAbsRefractory+...

                        round(rand*lengthAbsRefractory);

                    % disable spike potential for refractory period

                    % set all values in this range to 0

                    ANspikes(j,t+1:t+nPointsRefractory)=0;

                end

            end  %t

            % segment debugging

            % plotInstructions.figureNo=98;

            % plotInstructions.displaydt=ANdt;

            %  plotInstructions.numPlots=1;

            %  plotInstructions.subPlotNo=1;

            % UTIL_plotMatrix(ANspikes, plotInstructions);

            % and save it. NB, AN is now on 'speedUp' time

            ANoutput(:, reducedSegmentPTR: shorterSegmentEndPTR)=ANspikes;

            %% CN Macgregor first neucleus -------------------------------

            % input is from AN so ANdt is used throughout

            % Each CNneuron has a unique set of input fibers selected

            %  at random from the available AN fibers (CNinputfiberLists)

            % Create the dendritic current for that neuron

            % First get input spikes to this neuron

            synapseNo=1;

            for ch=1:nChannels

                for idx=1:nCNneuronsPerChannel

                    % determine candidate fibers for this unit

                    fibersUsed=CNinputfiberLists(synapseNo,:);

                    % ANpsth has a bin width of dt

                    %  (just a simple sum across fibers)

                    AN_PSTH(synapseNo,:) = ...

                        sum(ANspikes(fibersUsed,:), 1);

                    synapseNo=synapseNo+1;

                end

            end

            % One alpha function per spike

            [alphaRows alphaCols]=size(CNtrailingAlphas);

            for unitNo=1:nCNneurons

                CNcurrentTemp(unitNo,:)= ...

                    conv(AN_PSTH(unitNo,:),CNalphaFunction);

            end

            % 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);

            % 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-dts;

                    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))*dt;

                    dGk=-CN_Gk*dt./tauGk + CN_b*s;

                    dTh=-(CN_Th-CN_Th0)*dt/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)

                for t=1:reducedSegmentLength

                    CNtimeSinceLastSpike=CNtimeSinceLastSpike-dt;

                    s=CN_E>CN_Th0 & CNtimeSinceLastSpike<0 ;  % =1 if both conditions met

                    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))*dt;

                    dGk=-CN_Gk*dt./tauGk +CN_b*s;

                    CN_E=CN_E+dE;

                    CN_Gk=CN_Gk+dGk;

                    % add spike to CN_E and add resting potential (-60 mV)

                    CNmembranePotential(:,t)=CN_E+s.*(CN_Eb-CN_E)+CN_Er;

                end

            end

            % extract spikes.  A spike is a substantial upswing in voltage

            CN_spikes=CNmembranePotential> -0.01;

            % now remove any spike that is immediately followed by a spike

            % NB 'find' works on columns (whence the transposing)

            CN_spikes=CN_spikes';

            idx=find(CN_spikes);

            idx=idx(1:end-1);

            CN_spikes(idx+1)=0;

            CN_spikes=CN_spikes';

            % segment debugging

            % plotInstructions.figureNo=98;

            % plotInstructions.displaydt=ANdt;

            %  plotInstructions.numPlots=1;

            %  plotInstructions.subPlotNo=1;

            % UTIL_plotMatrix(CN_spikes, plotInstructions);

            % and save it

            CNoutput(:, reducedSegmentPTR:shorterSegmentEndPTR)=...

                CN_spikes;

            %% IC ----------------------------------------------

                %  MacGregor or some other second order neurons

                % combine CN neurons in same channel, i.e. same BF & same tauCa

                %  to generate inputs to single IC unit

                channelNo=0;

                for idx=1:nCNneuronsPerChannel:nCNneurons-nCNneuronsPerChannel+1;

                    channelNo=channelNo+1;

                    CN_PSTH(channelNo,:)=...

                        sum(CN_spikes(idx:idx+nCNneuronsPerChannel-1,:));

                end

                [alphaRows alphaCols]=size(ICtrailingAlphas);

                for ICneuronNo=1:nICcells

                    ICcurrentTemp(ICneuronNo,:)= ...

                        conv(CN_PSTH(ICneuronNo,:),  IC_CNalphaFunction);

                end

                % add the unused current from the previous convolution

                ICcurrentTemp(:,1:alphaCols)=ICcurrentTemp(:,1:alphaCols)...

                    + ICtrailingAlphas;

                % take what is required and keep the trailing part for next time

                inputCurrent=ICcurrentTemp(:, 1:reducedSegmentLength);

                ICtrailingAlphas=ICcurrentTemp(:, reducedSegmentLength+1:end);

                if IC_c==0

                    % faster computation when threshold is stable (C==0)

                    for t=1:reducedSegmentLength

                        s=IC_E>IC_Th0;

                        dE = (-IC_E/IC_tauM + inputCurrent(:,t)/IC_cap +...

                            (IC_Gk/IC_cap).*(IC_Ek-IC_E))*dt;

                        dGk=-IC_Gk*dt/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))*dt;

                        IC_E=IC_E+dE;

                        s=IC_E>IC_Th;

                        ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er;

                        dGk=-IC_Gk*dt/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)*dt/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

                    ICfiberTypeRates(tauCount, ...

                        reducedSegmentPTR:shorterSegmentEndPTR)=...

                        sum(ICspikes(firstCell:lastCell, :))...

                        /(nCellsPerTau*ANdt);

                    firstCell=firstCell+nCellsPerTau;

                    lastCell=lastCell+nCellsPerTau;

                end

                ICoutput(:, reducedSegmentPTR:shorterSegmentEndPTR)=ICspikes;

                if nBFs==1  % single channel

                    x= repmat(ICmembranePotential(1,:), ANspeedUpFactor,1);

                    x= reshape(x,1,segmentLength);

                    if nANfiberTypes>1  % save HSR and LSR

                        y= repmat(ICmembranePotential(end,:), ANspeedUpFactor,1);

                        y= reshape(y,1,segmentLength);

                        x=[x; y];

                    end

                    ICmembraneOutput(:, segmentStartPTR:segmentEndPTR)= x;

                end

                % estimate efferent effects.

                % ARis based on LSR units. LSR channels are 1:nBF

                if nANfiberTypes>1  % AR is multi-channel only

                    ARAttSeg=sum(ICspikes(1:nBFs,:),1)/ANdt;

                    [ARAttSeg, ARboundary] = ...

                        filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundary);

                    ARAttSeg=ARAttSeg-ARrateThreshold;

                    ARAttSeg(ARAttSeg<0)=0;   % prevent negative strengths

                    % scale up to dt from ANdt

                    x=    repmat(ARAttSeg, ANspeedUpFactor,1);

                    x=reshape(x,1,segmentLength);

                    ARattenuation(segmentStartPTR:segmentEndPTR)=...

                        (1-ARrateToAttenuationFactor* x);

                    ARattenuation(ARattenuation<0)=0.001;

                else

                    % single channel model; disable AR

                    ARattenuation(segmentStartPTR:segmentEndPTR)=...

                        ones(1,segmentLength);

                end

                % MOC attenuation using HSR response only

                % Separate MOC effect for each BF

                HSRbegins=nBFs*(nANfiberTypes-1)+1;

                rates=ICspikes(HSRbegins:end,:)/ANdt;

                for idx=1:nBFs

                    [smoothedRates, MOCboundary{idx}] = ...

                        filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ...

                        MOCboundary{idx});

                    MOCattSegment(idx,:)=smoothedRates;

                    % expand timescale back to model dt from ANdt

                    x= repmat(MOCattSegment(idx,:), ANspeedUpFactor,1);

                    x= reshape(x,1,segmentLength);

                    MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ...

                        (1- MOCrateToAttenuationFactor*  x);

                end

                MOCattenuation(MOCattenuation<0)=0.04;

                % segment debugging

                % plotInstructions.figureNo=98;

                % plotInstructions.displaydt=ANdt;

                %  plotInstructions.numPlots=1;

                %  plotInstructions.subPlotNo=1;

                % UTIL_plotMatrix(ICspikes, plotInstructions);

    end     % AN_spikesOrProbability

    segmentStartPTR=segmentStartPTR+segmentLength;

    reducedSegmentPTR=reducedSegmentPTR+reducedSegmentLength;

end  % segment

path(restorePath)