rmeddis@38: function test_binaural
rmeddis@38: % test_binaural is a first attempt to produce a binaural model
rmeddis@38: %  incorporating MSO and IC models.
rmeddis@38: % The monaural response is computed first for left and right stimuli
rmeddis@38: %  before using the CN response as input to the binaural MSO model
rmeddis@38: %  that, in turn, feeds a single cell IC model.
rmeddis@38: %
rmeddis@38: % The function has no arguments and everything is set up internally
rmeddis@38: %
rmeddis@38: % #1
rmeddis@38: % Identify the file (in 'MAPparamsName') containing the model parameters
rmeddis@38: %  the default is 'PL' which uses primary like neurons in the CN to
rmeddis@38: %  simulate spherical bushy cells
rmeddis@38: %
rmeddis@38: % #2
rmeddis@38: % Set AN_spikesOrProbability'). to 'spikes'
rmeddis@38: %
rmeddis@38: % #3
rmeddis@38: % Choose between a tone signal or file input (in 'signalType')
rmeddis@38: %
rmeddis@38: % #4
rmeddis@38: % Set the signal rms level (in leveldBSPL)
rmeddis@38: %
rmeddis@38: % #5
rmeddis@38: % Identify the channels in terms of their best frequencies in the vector
rmeddis@38: %  BFlist. This is currently a single-channel model, so only one BF needed
rmeddis@38: %
rmeddis@38: % #6
rmeddis@38: % Last minute changes to the parameters fetched earlier can be made using
rmeddis@38: %  the cell array of strings 'paramChanges'.
rmeddis@38: %  Each string must have the same format as the corresponding line in the
rmeddis@38: %  file identified in 'MAPparamsName'
rmeddis@38: % Currently this is used to specify that only HSR fibers are used and
rmeddis@38: %  for changing the current per AN spike at the CN dendrite
rmeddis@38: %
rmeddis@38: % #7
rmeddis@38: % specify the parameters of the MSO cells in the MSOParams structure
rmeddis@38: %
rmeddis@38: % #8
rmeddis@38: % specify the parameters of the IC cells in the ICMSOParams structure
rmeddis@38: %
rmeddis@38: % #9
rmeddis@38: % identify the plots required from MAP1_14 (i.e. before the bonaural model)
rmeddis@38: %
rmeddis@38: % #10
rmeddis@38: % Specify ITDs. The program cycles through different ITDs
rmeddis@38: %
rmeddis@38: 
rmeddis@38: global CNoutput dtSpikes
rmeddis@38: dbstop if error
rmeddis@38: restorePath=path;
rmeddis@38: addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'], ...
rmeddis@38:     ['..' filesep 'utilities'])
rmeddis@38: 
rmeddis@38: %%  #1 monaural model parameter file name
rmeddis@38: MAPparamsName='PL';
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% #2 'spikes' are mandatory for this model
rmeddis@38: AN_spikesOrProbability='spikes';
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% #3 pure tone, harmonic sequence or speech file input
rmeddis@38: signalType= 'tones';
rmeddis@38: sampleRate= 50000;
rmeddis@38: duration=0.050;                 % seconds
rmeddis@38: rampDuration=.005;              % raised cosine ramp (seconds)
rmeddis@38: beginSilence=0.050;
rmeddis@38: endSilence=0.050;
rmeddis@38: toneFrequency= 750;            % or a pure tone (Hz)
rmeddis@38: 
rmeddis@38: %   or
rmeddis@38: % harmonic sequence (Hz)
rmeddis@38: % F0=210;
rmeddis@38: % toneFrequency= F0:F0:8000;
rmeddis@38: 
rmeddis@38: %   or
rmeddis@38: % signalType= 'file';
rmeddis@38: % fileName='twister_44kHz';
rmeddis@38: 
rmeddis@38: if strcmp(signalType, 'file')
rmeddis@38:     % needed for labeling plot
rmeddis@38:     showMapOptions.fileName=fileName;
rmeddis@38: else
rmeddis@38:     showMapOptions.fileName=[];
rmeddis@38: end
rmeddis@38: 
rmeddis@38: %% #4 rms level
rmeddis@38: leveldBSPL= 70;
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% #5 number of channels in the model
rmeddis@38: BFlist=toneFrequency;
rmeddis@38: 
rmeddis@38: %% #6 change model parameters
rmeddis@38: 
rmeddis@38: paramChanges={'IHCpreSynapseParams.tauCa=80e-6;',...
rmeddis@38:     'MacGregorMultiParams.currentPerSpike=0.800e-6;'};
rmeddis@38: 
rmeddis@38: % Parameter changes can be used to change one or more model parameters
rmeddis@38: %  *after* the MAPparams file has been read
rmeddis@38: 
rmeddis@38: %% #7 MSO parameters
rmeddis@38: MSOParams.nNeuronsPerBF=	10;   % N neurons per BF
rmeddis@38: MSOParams.type = 'primary-like cell';
rmeddis@38: MSOParams.fibersPerNeuron=4;   % N input fibers
rmeddis@38: MSOParams.dendriteLPfreq=2000;  % dendritic filter
rmeddis@38: MSOParams.currentPerSpike=0.11e-6; % (A) per spike
rmeddis@38: MSOParams.currentPerSpike=0.5e-6; % (A) per spike
rmeddis@38: MSOParams.Cap=4.55e-9;   % cell capacitance (Siemens)
rmeddis@38: MSOParams.tauM=5e-4;     % membrane time constant (s)
rmeddis@38: MSOParams.Ek=-0.01;      % K+ eq. potential (V)
rmeddis@38: MSOParams.dGkSpike=3.64e-5; % K+ cond.shift on spike,S
rmeddis@38: MSOParams.tauGk=	0.0012; % K+ conductance tau (s)
rmeddis@38: MSOParams.Th0=	0.01;   % equilibrium threshold (V)
rmeddis@38: MSOParams.c=	0.01;       % threshold shift on spike, (V)
rmeddis@38: MSOParams.tauTh=	0.015;  % variable threshold tau
rmeddis@38: MSOParams.Er=-0.06;      % resting potential (V)
rmeddis@38: MSOParams.Eb=0.06;       % spike height (V)
rmeddis@38: MSOParams.debugging=0;        % (special)
rmeddis@38: MSOParams.wideband=0;         % special for wideband units
rmeddis@38: 
rmeddis@38: %% #8 IC parameters
rmeddis@38: ICchopperParams.nNeuronsPerBF=	10;   % N neurons per BF
rmeddis@38: ICchopperParams.type = 'chopper cell';
rmeddis@38: ICchopperParams.fibersPerNeuron=10;  % N input fibers
rmeddis@38: ICchopperParams.dendriteLPfreq=50;   % dendritic filter
rmeddis@38: ICchopperParams.currentPerSpike=50e-9; % *per spike
rmeddis@38: ICchopperParams.currentPerSpike=100e-9; % *per spike
rmeddis@38: ICchopperParams.Cap=1.67e-8; % ??cell capacitance (Siemens)
rmeddis@38: ICchopperParams.tauM=0.002;  % membrane time constant (s)
rmeddis@38: ICchopperParams.Ek=-0.01;    % K+ eq. potential (V)
rmeddis@38: ICchopperParams.dGkSpike=1.33e-4; % K+ cond.shift on spike,S
rmeddis@38: ICchopperParams.tauGk=	0.0005;% K+ conductance tau (s)
rmeddis@38: ICchopperParams.Th0=	0.01; % equilibrium threshold (V)
rmeddis@38: ICchopperParams.c=	0;        % threshold shift on spike, (V)
rmeddis@38: ICchopperParams.tauTh=	0.02; % variable threshold tau
rmeddis@38: ICchopperParams.Er=-0.06;    % resting potential (V)
rmeddis@38: ICchopperParams.Eb=0.06;     % spike height (V)
rmeddis@38: ICchopperParams.PSTHbinWidth=	1e-4;
rmeddis@38: 
rmeddis@38: %% #9 delare 'showMap' options to control graphical output
rmeddis@38: % this applies to the monaural input model only
rmeddis@38: showMapOptions.printModelParameters=0;   % prints all parameters
rmeddis@38: showMapOptions.showModelOutput=1;   % plot all stages if monaural input
rmeddis@38: 
rmeddis@38: %% #10 ITDs
rmeddis@38: % the program cycles through a range of stimulus ITDs
rmeddis@38: ITDs=0e-6:100e-6:2000e-6;
rmeddis@38: % ITDs=0; % single shot
rmeddis@38: 
rmeddis@38: %% Now start computing!
rmeddis@38: figure(98), clf, set(gcf, 'name', 'binaural demo')
rmeddis@38: MSOcounts=zeros(1,length(ITDs));
rmeddis@38: ICcounts=zeros(1,length(ITDs));
rmeddis@38: ITDcount=0;
rmeddis@38: for ITD=ITDs
rmeddis@38:     ITDcount=ITDcount+1;
rmeddis@38:     delaySamples=round(ITD* sampleRate);
rmeddis@38:     %% Generate stimuli
rmeddis@38:     switch signalType
rmeddis@38:         case 'tones'
rmeddis@38:             % Create pure tone stimulus
rmeddis@38:             dt=1/sampleRate; % seconds
rmeddis@38:             time=dt: dt: duration;
rmeddis@38:             inputSignal=sum(sin(2*pi*toneFrequency'*time), 1);
rmeddis@38:             amp=10^(leveldBSPL/20)*28e-6;   % converts to Pascals (peak)
rmeddis@38:             inputSignal=amp*inputSignal;
rmeddis@38:             % apply ramps
rmeddis@38:             % catch rampTime error
rmeddis@38:             if rampDuration>0.5*duration, rampDuration=duration/2; end
rmeddis@38:             rampTime=dt:dt:rampDuration;
rmeddis@38:             ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
rmeddis@38:                 ones(1,length(time)-length(rampTime))];
rmeddis@38:             inputSignal=inputSignal.*ramp;
rmeddis@38:             ramp=fliplr(ramp);
rmeddis@38:             inputSignal=inputSignal.*ramp;
rmeddis@38:             % add silence
rmeddis@38:             intialSilence= zeros(1,round(beginSilence/dt));
rmeddis@38:             finalSilence= zeros(1,round(endSilence/dt));
rmeddis@38:             inputSignal= [intialSilence inputSignal finalSilence];
rmeddis@38: 
rmeddis@38:         case 'file'
rmeddis@38:             %% file input simple or mixed
rmeddis@38:             [inputSignal sampleRate]=wavread(fileName);
rmeddis@38:             dt=1/sampleRate;
rmeddis@38:             inputSignal=inputSignal(:,1);
rmeddis@38:             targetRMS=20e-6*10^(leveldBSPL/20);
rmeddis@38:             rms=(mean(inputSignal.^2))^0.5;
rmeddis@38:             amp=targetRMS/rms;
rmeddis@38:             inputSignal=inputSignal*amp;
rmeddis@38:             intialSilence= zeros(1,round(0.1/dt));
rmeddis@38:             finalSilence= zeros(1,round(0.2/dt));
rmeddis@38:             inputSignal= [intialSilence inputSignal' finalSilence];
rmeddis@38:     end
rmeddis@38: 
rmeddis@38:     %% run the monaural model twice
rmeddis@38:     t=dt:dt:dt*length(inputSignal);
rmeddis@38:     for ear={'left','right'}
rmeddis@38:         figure(98), subplot(4,1,1), colour='b'; hold off
rmeddis@38:         switch ear{1}
rmeddis@38:             case 'right'
rmeddis@38:                 inputSignal=circshift(inputSignal', delaySamples)';
rmeddis@38:                 hold on
rmeddis@38:                 colour='r';
rmeddis@38:         end
rmeddis@38:         plot(t, inputSignal, colour)
rmeddis@38:         title ('binaural inputs signals')
rmeddis@38:         ylabel('Pa'), xlabel('time')
rmeddis@38:         xlim([0 max(t)])
rmeddis@38: 
rmeddis@38:         % call to monaural model
rmeddis@38:         MAP1_14(inputSignal, sampleRate, BFlist, ...
rmeddis@38:             MAPparamsName, AN_spikesOrProbability, paramChanges);
rmeddis@38: 
rmeddis@38:         % the model run is now complete. Now display the results
rmeddis@38:         UTIL_showMAP(showMapOptions, paramChanges)
rmeddis@38: 
rmeddis@38:         % copy the CN inputspiking pattern to the binaural display
rmeddis@38:         figure(98)
rmeddis@38:         switch ear{1}
rmeddis@38:             case 'left'
rmeddis@38:                 CNoutputLeft=CNoutput;
rmeddis@38:                 subplot(4,2,3)
rmeddis@38:             case 'right'
rmeddis@38:                 CNoutputRight=CNoutput;
rmeddis@38:                 subplot(4,2,4)
rmeddis@38:         end
rmeddis@38:         plotInstructions=[];
rmeddis@38:         plotInstructions.axes=gca;
rmeddis@38:         plotInstructions.displaydt=dtSpikes;
rmeddis@38:         plotInstructions.title= ['CN spikes ' ear{1}];
rmeddis@38:         plotInstructions.rasterDotSize=2;
rmeddis@38:         if sum(sum(CNoutput))<100
rmeddis@38:             plotInstructions.rasterDotSize=3;
rmeddis@38:         end
rmeddis@38:         UTIL_plotMatrix(CNoutput, plotInstructions);
rmeddis@38: 
rmeddis@38:     end % left/ right ear
rmeddis@38: 
rmeddis@38:     %%  MSO
rmeddis@38:     % run MSO model using left and right CN spikes
rmeddis@38:     MSOspikes=MSO(CNoutputLeft,CNoutputRight, dtSpikes, MSOParams);
rmeddis@38:     
rmeddis@38:     sumspikes=sum(sum(MSOspikes));
rmeddis@38:     disp(['ITD/ MSO spikes count= ' num2str([ITD sumspikes])])
rmeddis@38:     MSOcounts(ITDcount)=sumspikes;
rmeddis@38:     figure(98), subplot(4,2, 8), cla, hold off, plot(ITDs*1e6, MSOcounts)
rmeddis@38: 
rmeddis@38:     %% IC chopper
rmeddis@38:     % run IC model using all MSO spikes as input to a single IC cell
rmeddis@38:     ICMSOspikes=ICchopper(MSOspikes, dtSpikes, ICchopperParams);
rmeddis@38:     
rmeddis@38:     sumspikes=sum(sum(ICMSOspikes));
rmeddis@38:     disp(['ITD/ ICMSO spikes count= ' num2str([ITD sumspikes])])
rmeddis@38:     ICcounts(ITDcount)=sumspikes;
rmeddis@38:     figure(98), subplot(4,2,8),hold on, plot(ITDs*1e6, ICcounts,'r')
rmeddis@38:     xlabel('ITD'), ylabel(' spike count')
rmeddis@38:     title('MSO (blue)/ IC (red) spike counts')
rmeddis@38:     legend({'MSO','IC'})
rmeddis@38: 
rmeddis@38: end % ITDs
rmeddis@38: 
rmeddis@38: path(restorePath)
rmeddis@38: 
rmeddis@38: 
rmeddis@38: function MSOspikes=MSO(CNoutputLeft,CNoutputRight, dtSpikes, MSOparams)
rmeddis@38: %%
rmeddis@38: [nMSOcells nEpochs]=size(CNoutputLeft);
rmeddis@38: inputCurrent=zeros(nMSOcells, nEpochs);
rmeddis@38: MSOmembranePotential=inputCurrent;
rmeddis@38: 
rmeddis@38: MSO_tauM=MSOparams.tauM;
rmeddis@38: MSO_tauGk=MSOparams.tauGk;
rmeddis@38: MSO_tauTh=MSOparams.tauTh;
rmeddis@38: MSO_cap=MSOparams.Cap;
rmeddis@38: MSO_c=MSOparams.c;
rmeddis@38: MSO_b=MSOparams.dGkSpike;
rmeddis@38: MSO_Th0=MSOparams.Th0;
rmeddis@38: MSO_Ek=MSOparams.Ek;
rmeddis@38: MSO_Eb= MSOparams.Eb;
rmeddis@38: MSO_Er=MSOparams.Er;
rmeddis@38: 
rmeddis@38: MSO_E=zeros(nMSOcells,1);
rmeddis@38: MSO_Gk=zeros(nMSOcells,1);
rmeddis@38: MSO_Th=MSO_Th0*ones(nMSOcells,1);
rmeddis@38: 
rmeddis@38: % Dendritic filtering, all spikes are replaced by CNalphaFunction functions
rmeddis@38: MSOdendriteLPfreq= MSOparams.dendriteLPfreq;
rmeddis@38: MSOcurrentPerSpike=MSOparams.currentPerSpike;
rmeddis@38: MSOspikeToCurrentTau=1/(2*pi*MSOdendriteLPfreq);
rmeddis@38: t=dtSpikes:dtSpikes:5*MSOspikeToCurrentTau;
rmeddis@38: MSO_CNalphaFunction= (MSOcurrentPerSpike / ...
rmeddis@38:     MSOspikeToCurrentTau)*t.*exp(-t / MSOspikeToCurrentTau);
rmeddis@38: % show alpha function
rmeddis@38: % figure(84), subplot(4,2,2), plot(t,MSO_CNalphaFunction)
rmeddis@38: % title(['LP cutoff ' num2str(MSOdendriteLPfreq)])
rmeddis@38: 
rmeddis@38: % convert CN spikes to post-dendritic current
rmeddis@38: CN_spikes=CNoutputLeft+CNoutputRight;
rmeddis@38: for i=1:nMSOcells
rmeddis@38:     x= conv2(CN_spikes(i,:),  MSO_CNalphaFunction);
rmeddis@38:     inputCurrent(i,:)=x(1:nEpochs);
rmeddis@38: end
rmeddis@38: 
rmeddis@38: if MSO_c==0
rmeddis@38:     % faster computation when threshold is stable (c==0)
rmeddis@38:     for t=1:nEpochs
rmeddis@38:         s=MSO_E>MSO_Th0;
rmeddis@38:         dE = (-MSO_E/MSO_tauM + inputCurrent(:,t)/MSO_cap +...
rmeddis@38:             (MSO_Gk/MSO_cap).*(MSO_Ek-MSO_E))*dtSpikes;
rmeddis@38:         dGk=-MSO_Gk*dtSpikes/MSO_tauGk +MSO_b*s;
rmeddis@38:         MSO_E=MSO_E+dE;
rmeddis@38:         MSO_Gk=MSO_Gk+dGk;
rmeddis@38:         MSOmembranePotential(:,t)=MSO_E+s.*(MSO_Eb-MSO_E)+MSO_Er;
rmeddis@38:     end
rmeddis@38: else
rmeddis@38:     %  threshold is changing (MSO_c>0; e.g. bushy cell)
rmeddis@38:     for t=1:nEpochs
rmeddis@38:         dE = (-MSO_E/MSO_tauM + ...
rmeddis@38:             inputCurrent(:,t)/MSO_cap + (MSO_Gk/MSO_cap)...
rmeddis@38:             .*(MSO_Ek-MSO_E))*dtSpikes;
rmeddis@38:         MSO_E=MSO_E+dE;
rmeddis@38:         s=MSO_E>MSO_Th;
rmeddis@38:         MSOmembranePotential(:,t)=MSO_E+s.*(MSO_Eb-MSO_E)+MSO_Er;
rmeddis@38:         dGk=-MSO_Gk*dtSpikes/MSO_tauGk +MSO_b*s;
rmeddis@38:         MSO_Gk=MSO_Gk+dGk;
rmeddis@38: 
rmeddis@38:         % After a spike, the threshold is raised
rmeddis@38:         % otherwise it settles to its baseline
rmeddis@38:         dTh=-(MSO_Th-MSO_Th0)*dtSpikes/MSO_tauTh +s*MSO_c;
rmeddis@38:         MSO_Th=MSO_Th+dTh;
rmeddis@38:     end
rmeddis@38: end
rmeddis@38: 
rmeddis@38: figure(98),subplot(4,1,3)
rmeddis@38: hold off, imagesc(MSOmembranePotential)
rmeddis@38: title ('MSO (V)')
rmeddis@38: 
rmeddis@38: MSOspikes=MSOmembranePotential> -0.01;
rmeddis@38: % Remove any spike that is immediately followed by a spike
rmeddis@38: %  NB 'find' works on columns (whence the transposing)
rmeddis@38: MSOspikes=MSOspikes';
rmeddis@38: idx=find(MSOspikes);
rmeddis@38: idx=idx(1:end-1);
rmeddis@38: MSOspikes(idx+1)=0;
rmeddis@38: MSOspikes=MSOspikes';
rmeddis@38: 
rmeddis@38: 
rmeddis@38: function ICMSOspikes=ICchopper(ICMSOspikes, dtSpikes, ICMSOParams)
rmeddis@38: %%
rmeddis@38: ICMSOspikes=sum(ICMSOspikes);
rmeddis@38: [nICMSOcells nEpochs]=size(ICMSOspikes);
rmeddis@38: inputCurrent=zeros(nICMSOcells, nEpochs);
rmeddis@38: ICMSOmembranePotential=inputCurrent;
rmeddis@38: 
rmeddis@38: ICMSO_tauM=ICMSOParams.tauM;
rmeddis@38: ICMSO_tauGk=ICMSOParams.tauGk;
rmeddis@38: ICMSO_tauTh=ICMSOParams.tauTh;
rmeddis@38: ICMSO_cap=ICMSOParams.Cap;
rmeddis@38: ICMSO_c=ICMSOParams.c;
rmeddis@38: ICMSO_b=ICMSOParams.dGkSpike;
rmeddis@38: ICMSO_Th0=ICMSOParams.Th0;
rmeddis@38: ICMSO_Ek=ICMSOParams.Ek;
rmeddis@38: ICMSO_Eb= ICMSOParams.Eb;
rmeddis@38: ICMSO_Er=ICMSOParams.Er;
rmeddis@38: 
rmeddis@38: ICMSO_E=zeros(nICMSOcells,1);
rmeddis@38: ICMSO_Gk=zeros(nICMSOcells,1);
rmeddis@38: ICMSO_Th=ICMSO_Th0*ones(nICMSOcells,1);
rmeddis@38: 
rmeddis@38: % Dendritic filtering, all spikes are replaced by CNalphaFunction functions
rmeddis@38: ICMSOdendriteLPfreq= ICMSOParams.dendriteLPfreq;
rmeddis@38: ICMSOcurrentPerSpike=ICMSOParams.currentPerSpike;
rmeddis@38: ICMSOspikeToCurrentTau=1/(2*pi*ICMSOdendriteLPfreq);
rmeddis@38: t=dtSpikes:dtSpikes:5*ICMSOspikeToCurrentTau;
rmeddis@38: ICMSOalphaFunction= (ICMSOcurrentPerSpike / ...
rmeddis@38:     ICMSOspikeToCurrentTau)*t.*exp(-t / ICMSOspikeToCurrentTau);
rmeddis@38: % show alpha function
rmeddis@38: % figure(84), subplot(4,2,5), plot(t,ICMSOalphaFunction)
rmeddis@38: % title(['IC MSO LP cutoff ' num2str(ICMSOdendriteLPfreq)])
rmeddis@38: 
rmeddis@38: % post-dendritic current
rmeddis@38: for i=1:nICMSOcells
rmeddis@38:     x= conv2(1*ICMSOspikes(i,:),  ICMSOalphaFunction);
rmeddis@38:     inputCurrent(i,:)=x(1:nEpochs);
rmeddis@38: end
rmeddis@38: 
rmeddis@38: if ICMSO_c==0
rmeddis@38:     % faster computation when threshold is stable (c==0)
rmeddis@38:     for t=1:nEpochs
rmeddis@38:         s=ICMSO_E>ICMSO_Th0;
rmeddis@38:         dE = (-ICMSO_E/ICMSO_tauM + inputCurrent(:,t)/ICMSO_cap +...
rmeddis@38:             (ICMSO_Gk/ICMSO_cap).*(ICMSO_Ek-ICMSO_E))*dtSpikes;
rmeddis@38:         dGk=-ICMSO_Gk*dtSpikes/ICMSO_tauGk +ICMSO_b*s;
rmeddis@38:         ICMSO_E=ICMSO_E+dE;
rmeddis@38:         ICMSO_Gk=ICMSO_Gk+dGk;
rmeddis@38:         ICMSOmembranePotential(:,t)=ICMSO_E+s.*(ICMSO_Eb-ICMSO_E)+ICMSO_Er;
rmeddis@38:     end
rmeddis@38: else
rmeddis@38:     %  threshold is changing (ICMSO_c>0; e.g. bushy cell)
rmeddis@38:     for t=1:nEpochs
rmeddis@38:         dE = (-ICMSO_E/ICMSO_tauM + ...
rmeddis@38:             inputCurrent(:,t)/ICMSO_cap + (ICMSO_Gk/ICMSO_cap)...
rmeddis@38:             .*(ICMSO_Ek-ICMSO_E))*dtSpikes;
rmeddis@38:         ICMSO_E=ICMSO_E+dE;
rmeddis@38:         s=ICMSO_E>ICMSO_Th;
rmeddis@38:         ICMSOmembranePotential(:,t)=ICMSO_E+s.*(ICMSO_Eb-ICMSO_E)+ICMSO_Er;
rmeddis@38:         dGk=-ICMSO_Gk*dtSpikes/ICMSO_tauGk +ICMSO_b*s;
rmeddis@38:         ICMSO_Gk=ICMSO_Gk+dGk;
rmeddis@38: 
rmeddis@38:         % After a spike, the threshold is raised
rmeddis@38:         % otherwise it settles to its baseline
rmeddis@38:         dTh=-(ICMSO_Th-ICMSO_Th0)*dtSpikes/ICMSO_tauTh +s*ICMSO_c;
rmeddis@38:         ICMSO_Th=ICMSO_Th+dTh;
rmeddis@38:     end
rmeddis@38: end
rmeddis@38: 
rmeddis@38: t=dtSpikes:dtSpikes:dtSpikes*length(ICMSOmembranePotential);
rmeddis@38: figure(98),subplot(4,2,7)
rmeddis@38: plot(t, ICMSOmembranePotential)
rmeddis@38: ylim([-0.07 0]), xlim([0 max(t)])
rmeddis@38: title('IC (V)')
rmeddis@38: 
rmeddis@38: ICMSOspikes=ICMSOmembranePotential> -0.01;
rmeddis@38: % now remove any spike that is immediately followed by a spike
rmeddis@38: % NB 'find' works on columns (whence the transposing)
rmeddis@38: ICMSOspikes=ICMSOspikes';
rmeddis@38: idx=find(ICMSOspikes);
rmeddis@38: idx=idx(1:end-1);
rmeddis@38: ICMSOspikes(idx+1)=0;
rmeddis@38: ICMSOspikes=ICMSOspikes';
rmeddis@38: