MAP

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  savedParamChanges savedBFlist saveAN_spikesOrProbability ...

    saveMAPparamsName savedInputSignal dt dtSpikes ...

    OMEextEarPressure TMoutput OMEoutput DRNLoutput...

    IHC_cilia_output IHCrestingCiliaCond IHCrestingV...

    IHCoutput ANprobRateOutput ANoutput savePavailable saveNavailable ...

    ANtauCas CNtauGk CNoutput ICoutput ICmembraneOutput ICfiberTypeRates...

    MOCattenuation ARattenuation

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

% needed by the user; so the following will suffice

%   global dtSpikes 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 dtSpikes).

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

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

savedParamChanges=paramChanges;

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=ceil(sampleRate/AN_IHCsynapseParams.spikesTargetSampleRate);

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

% spiking activity is computed at longer sampling intervals (dtSpikes)

%  so it has a smaller number of epochs per segment(see 'ANspeeUpFactor' above)

% AN CN and IC all use this sample interval

dtSpikes=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.OMEstapesHPcutoff;

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;

minMOCattenuation=10^(DRNLParams.minMOCattenuationdB/20);

% smoothing filter for MOC

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

b0=1+ a1;

MOCfilt_b=b0; MOCfilt_a=[a0 a1];

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

b0=1+ a1;

MOCfiltProb_b=b0; MOCfiltProb_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;

% DRNLa2=DRNLParams.a2;

% DRNLb=DRNLParams.b;

DRNLc=DRNLParams.c;

linGAIN=DRNLParams.g;

ctBM=10e-9*10^(DRNLParams.ctBMdB/20);

CtS=ctBM/DRNLa;

%

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

% i.e. b= [1  dt/tc-1]  and a= dt/IHC_cilia_RPParams.tc

% 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 vector (one per 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);

% only one stream of available transmitter will be saved

saveNavailableSeg=zeros(1,reducedSegmentLength);

saveNavailable=zeros(1,reducedSignalLength);

% spikes     % !  !  !    ! !        !   !  !

lengthAbsRefractory= round(AN_refractory_period/dtSpikes);

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

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

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

AN_rdt= repmat(AN_IHCsynapseParams.r*dtSpikes, 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=dtSpikes:dtSpikes:5*CNspikeToCurrentTau;

CNalphaFunction= (1 / ...

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

CNalphaFunction=CNalphaFunction*CNcurrentPerSpike;

% figure(98), plot(t,CNalphaFunction), xlim([0 .020]), xlabel('time (s)'), ylabel('I')

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

% 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=dtSpikes:dtSpikes: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

%         figure(88), plot(MOC)

        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

        % Nick's compression algorithm

        abs_x= abs(nonlinOutput);

        signs= sign(nonlinOutput);

        belowThreshold= abs_x<CtS;

        nonlinOutput(belowThreshold)= DRNLa *nonlinOutput(belowThreshold);

        aboveThreshold=~belowThreshold;

        nonlinOutput(aboveThreshold)= signs(aboveThreshold) *ctBM .* ...

            exp(DRNLc *log( DRNLa*abs_x(aboveThreshold)/ctBM ));

%         %    original   broken stick instantaneous compression

%         holdY=nonlinOutput;

%         abs_x = abs(nonlinOutput);

%         nonlinOutput=sign(nonlinOutput).*min(DRNLa*abs_x, DRNLb*abs_x.^DRNLc);

%

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

%         end

%         nonlinOutput=y;

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

    %     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 to obtain the driving voltage

    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 channels open - apply time membrane constant

    for idx=1:segmentLength

        % mICaINF is the current 'target' value of mICa

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

        mICa(:,idx)=mICaCurrent;

    end

    % calcium current

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

    % apply calcium channel time constant

    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'

            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) :dt: tnow

               %    pFired= 1 - II(1-p(t)),

               % we need a running account of cumProb=II(1-p(t)) in order

               % not to have to recompute this for each value of 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 acoustic reflex efferent effect:  0 < ARattenuation > 1

            [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(MOCfiltProb_b, MOCfiltProb_a, rates(idx,:), ...

                        MOCprobBoundary{idx});

                    smoothedRates=smoothedRates-MOCrateThresholdProb;

                    smoothedRates=max(smoothedRates, 0);

                    x=(1- smoothedRates* rateToAttenuationFactorProb);

                    MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= x;

                end

            end

            MOCattenuation(MOCattenuation<minMOCattenuation)= minMOCattenuation;

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

                % 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

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

                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;

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

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

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

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

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

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

                AN_available = AN_available + replenish - ejected ...

                    + reprocessed;

                saveNavailableSeg(:,ANtimeCount)=AN_available(end,:); % only last channel

                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=dtSpikes;

            %  plotInstructions.numPlots=1;

            %  plotInstructions.subPlotNo=1;

            % UTIL_plotMatrix(ANspikes, plotInstructions);

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

            ANoutput(:, reducedSegmentPTR: shorterSegmentEndPTR)=ANspikes;

            % monitor synapse contents (only sometimes used)

            saveNavailable(reducedSegmentPTR:reducedSegmentPTR+reducedSegmentLength-1)=...

                saveNavailableSeg;

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

            % input is from AN so dtSpikes 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 dtSpikes

                    %  (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-dtSpikes;

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

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

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

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

                    dGk=-CN_Gk*dtSpikes./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=dtSpikes;

            %  plotInstructions.numPlots=1;

            %  plotInstructions.subPlotNo=1;