Mercurial > hg > map

comparison userProgramsTim/MAP1_14_olde_version_dont_use.m @ 38:c2204b18f4a2 tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
End nov big change
author Ray Meddis <rmeddis@essex.ac.uk>
date Mon, 28 Nov 2011 13:34:28 +0000
parents
children
comparison
equal deleted inserted replaced
37:771a643d5c29 38:c2204b18f4a2
1 function MAP1_14(inputSignal, sampleRate, BFlist, MAPparamsName, ...
2 AN_spikesOrProbability, paramChanges)
3 % To test this function use test_MAP1_14 in this folder
4 %
5 % example:
6 % <navigate to 'MAP1_14\MAP'>
7 % [inputSignal FS] = wavread('../wavFileStore/twister_44kHz');
8 % MAP1_14(inputSignal, FS, -1, 'Normal', 'probability', [])
9 %
10 % All arguments are mandatory.
11 %
12 % BFlist is a vector of BFs but can be '-1' to allow MAPparams to choose
13 % MAPparamsName='Normal'; % source of model parameters
14 % AN_spikesOrProbability='spikes'; % or 'probability'
15 % paramChanges is a cell array of strings that can be used to make last
16 % minute parameter changes, e.g., to simulate OHC loss
17 % e.g. paramChanges{1}= 'DRNLParams.a=0;'; % disable OHCs
18 % e.g. paramchanges={}; % no changes
19 % The model parameters are established in the MAPparams<***> file
20 % and stored as global
21
22 restorePath=path;
23 addpath (['..' filesep 'parameterStore'])
24
25 global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams
26 global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams
27
28 % All of the results of this function are stored as global
29 global dt ANdt savedBFlist saveAN_spikesOrProbability saveMAPparamsName...
30 savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ...
31 DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
32 IHCoutput ANprobRateOutput ANoutput savePavailable ANtauCas ...
33 CNtauGk CNoutput ICoutput ICmembraneOutput ICfiberTypeRates ...
34 MOCattenuation
35
36 % Normally only ICoutput(logical spike matrix) or ANprobRateOutput will be
37 % needed by the user; so the following will suffice
38 % global ANdt ICoutput ANprobRateOutput
39
40 % Note that sampleRate has not changed from the original function call and
41 % ANprobRateOutput is sampled at this rate
42 % However ANoutput, CNoutput and IC output are stored as logical
43 % 'spike' matrices using a lower sample rate (see ANdt).
44
45 % When AN_spikesOrProbability is set to probability,
46 % no spike matrices are computed.
47 % When AN_spikesOrProbability is set to 'spikes',
48 % no probability output is computed
49
50 % Efferent control variables are ARattenuation and MOCattenuation
51 % These are scalars between 1 (no attenuation) and 0.
52 % They are represented with dt=1/sampleRate (not ANdt)
53 % They are computed using either AN probability rate output
54 % or IC (spikes) output as approrpriate.
55 % AR is computed using across channel activity
56 % MOC is computed on a within-channel basis.
57
58 if nargin<1
59 error(' MAP1_14 is not a script but a function that must be called')
60 end
61
62 if nargin<6
63 paramChanges=[];
64 end
65 % Read parameters from MAPparams<***> file in 'parameterStore' folder
66 % Beware, 'BFlist=-1' is a legitimate argument for MAPparams<>
67 % It means that the calling program allows MAPparams to specify the list
68 cmd=['method=MAPparams' MAPparamsName ...
69 '(BFlist, sampleRate, 0, paramChanges);'];
70 eval(cmd);
71 BFlist=DRNLParams.nonlinCFs;
72
73 % save as global for later plotting if required
74 savedBFlist=BFlist;
75 saveAN_spikesOrProbability=AN_spikesOrProbability;
76 saveMAPparamsName=MAPparamsName;
77
78 dt=1/sampleRate;
79 duration=length(inputSignal)/sampleRate;
80 % segmentDuration is specified in parameter file (must be >efferent delay)
81 segmentDuration=method.segmentDuration;
82 segmentLength=round(segmentDuration/ dt);
83 segmentTime=dt*(1:segmentLength); % used in debugging plots
84
85 % all spiking activity is computed using longer epochs
86 ANspeedUpFactor=AN_IHCsynapseParams.ANspeedUpFactor; % e.g.5 times
87
88 % inputSignal must be row vector
89 [r c]=size(inputSignal);
90 if r>c, inputSignal=inputSignal'; end % transpose
91 % ignore stereo signals
92 inputSignal=inputSignal(1,:); % drop any second channel
93 savedInputSignal=inputSignal;
94
95 % Segment the signal
96 % The sgment length is given but the signal length must be adjusted to be a
97 % multiple of both the segment length and the reduced segmentlength
98 [nSignalRows signalLength]=size(inputSignal);
99 segmentLength=ceil(segmentLength/ANspeedUpFactor)*ANspeedUpFactor;
100 % Make the signal length a whole multiple of the segment length
101 nSignalSegments=ceil(signalLength/segmentLength);
102 padSize=nSignalSegments*segmentLength-signalLength;
103 pad=zeros(nSignalRows,padSize);
104 inputSignal=[inputSignal pad];
105 [ignore signalLength]=size(inputSignal);
106
107 % AN (spikes) is computed at a lower sample rate when spikes required
108 % so it has a reduced segment length (see 'ANspeeUpFactor' above)
109 % AN CN and IC all use this sample interval
110 ANdt=dt*ANspeedUpFactor;
111 reducedSegmentLength=round(segmentLength/ANspeedUpFactor);
112 reducedSignalLength= round(signalLength/ANspeedUpFactor);
113
114 %% Initialise with respect to each stage before computing
115 % by allocating memory,
116 % by computing constants
117 % by establishing easy to read variable names
118 % The computations are made in segments and boundary conditions must
119 % be established and stored. These are found in variables with
120 % 'boundary' or 'bndry' in the name
121
122 %% OME ---
123 % external ear resonances
124 OMEexternalResonanceFilters=OMEParams.externalResonanceFilters;
125 [nOMEExtFilters c]=size(OMEexternalResonanceFilters);
126 % details of external (outer ear) resonances
127 OMEgaindBs=OMEexternalResonanceFilters(:,1);
128 OMEgainScalars=10.^(OMEgaindBs/20);
129 OMEfilterOrder=OMEexternalResonanceFilters(:,2);
130 OMElowerCutOff=OMEexternalResonanceFilters(:,3);
131 OMEupperCutOff=OMEexternalResonanceFilters(:,4);
132 % external resonance coefficients
133 ExtFilter_b=cell(nOMEExtFilters,1);
134 ExtFilter_a=cell(nOMEExtFilters,1);
135 for idx=1:nOMEExtFilters
136 Nyquist=sampleRate/2;
137 [b, a] = butter(OMEfilterOrder(idx), ...
138 [OMElowerCutOff(idx) OMEupperCutOff(idx)]...
139 /Nyquist);
140 ExtFilter_b{idx}=b;
141 ExtFilter_a{idx}=a;
142 end
143 OMEExtFilterBndry=cell(2,1);
144 OMEextEarPressure=zeros(1,signalLength); % pressure at tympanic membrane
145
146 % pressure to velocity conversion using smoothing filter (50 Hz cutoff)
147 tau=1/(2*pi*50);
148 a1=dt/tau-1; a0=1;
149 b0=1+ a1;
150 TMdisp_b=b0; TMdisp_a=[a0 a1];
151 % figure(9), freqz(TMdisp_b, TMdisp_a)
152 OME_TMdisplacementBndry=[];
153
154 % OME high pass (simulates poor low frequency stapes response)
155 OMEhighPassHighCutOff=OMEParams.OMEstapesLPcutoff;
156 Nyquist=sampleRate/2;
157 [stapesDisp_b,stapesDisp_a] = butter(1, OMEhighPassHighCutOff/Nyquist, 'high');
158 % figure(10), freqz(stapesDisp_b, stapesDisp_a)
159
160 OMEhighPassBndry=[];
161
162 % OMEampStapes might be reducdant (use OMEParams.stapesScalar)
163 stapesScalar= OMEParams.stapesScalar;
164
165 % Acoustic reflex
166 efferentDelayPts=round(OMEParams.ARdelay/dt);
167 % smoothing filter
168 a1=dt/OMEParams.ARtau-1; a0=1;
169 b0=1+ a1;
170 ARfilt_b=b0; ARfilt_a=[a0 a1];
171
172 ARattenuation=ones(1,signalLength);
173 ARrateThreshold=OMEParams.ARrateThreshold; % may not be used
174 ARrateToAttenuationFactor=OMEParams.rateToAttenuationFactor;
175 ARrateToAttenuationFactorProb=OMEParams.rateToAttenuationFactorProb;
176 ARboundary=[];
177 ARboundaryProb=0;
178
179 % save complete OME record (stapes displacement)
180 OMEoutput=zeros(1,signalLength);
181 TMoutput=zeros(1,signalLength);
182
183 %% BM ---
184 % BM is represented as a list of locations identified by BF
185 DRNL_BFs=BFlist;
186 nBFs= length(DRNL_BFs);
187
188 % DRNLchannelParameters=DRNLParams.channelParameters;
189 DRNLresponse= zeros(nBFs, segmentLength);
190
191 MOCrateToAttenuationFactor=DRNLParams.rateToAttenuationFactor;
192 rateToAttenuationFactorProb=DRNLParams.rateToAttenuationFactorProb;
193 MOCrateThresholdProb=DRNLParams.MOCrateThresholdProb;
194
195 % smoothing filter for MOC
196 a1=dt/DRNLParams.MOCtau-1; a0=1;
197 b0=1+ a1;
198 MOCfilt_b=b0; MOCfilt_a=[a0 a1];
199 % figure(9), freqz(stapesDisp_b, stapesDisp_a)
200 MOCboundary=cell(nBFs,1);
201 MOCprobBoundary=cell(nBFs,1);
202
203 MOCattSegment=zeros(nBFs,reducedSegmentLength);
204 MOCattenuation=ones(nBFs,signalLength);
205
206 % if DRNLParams.a>0
207 % DRNLcompressionThreshold=10^((1/(1-DRNLParams.c))* ...
208 % log10(DRNLParams.b/DRNLParams.a));
209 % else
210 % DRNLcompressionThreshold=inf;
211 % end
212 DRNLcompressionThreshold=DRNLParams.cTh;
213 DRNLlinearOrder= DRNLParams.linOrder;
214 DRNLnonlinearOrder= DRNLParams.nonlinOrder;
215
216 DRNLa=DRNLParams.a;
217 DRNLb=DRNLParams.b;
218 DRNLc=DRNLParams.c;
219 linGAIN=DRNLParams.g;
220 %
221 % gammatone filter coefficients for linear pathway
222 bw=DRNLParams.linBWs';
223 phi = 2 * pi * bw * dt;
224 cf=DRNLParams.linCFs';
225 theta = 2 * pi * cf * dt;
226 cos_theta = cos(theta);
227 sin_theta = sin(theta);
228 alpha = -exp(-phi).* cos_theta;
229 b0 = ones(nBFs,1);
230 b1 = 2 * alpha;
231 b2 = exp(-2 * phi);
232 z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i;
233 z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i;
234 z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i;
235 tf = (z2 + z3) ./ z1;
236 a0 = abs(tf);
237 a1 = alpha .* a0;
238 GTlin_a = [b0, b1, b2];
239 GTlin_b = [a0, a1];
240 GTlinOrder=DRNLlinearOrder;
241 GTlinBdry=cell(nBFs,GTlinOrder);
242
243 % nonlinear gammatone filter coefficients
244 bw=DRNLParams.nlBWs';
245 phi = 2 * pi * bw * dt;
246 cf=DRNLParams.nonlinCFs';
247 theta = 2 * pi * cf * dt;
248 cos_theta = cos(theta);
249 sin_theta = sin(theta);
250 alpha = -exp(-phi).* cos_theta;
251 b0 = ones(nBFs,1);
252 b1 = 2 * alpha;
253 b2 = exp(-2 * phi);
254 z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i;
255 z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i;
256 z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i;
257 tf = (z2 + z3) ./ z1;
258 a0 = abs(tf);
259 a1 = alpha .* a0;
260 GTnonlin_a = [b0, b1, b2];
261 GTnonlin_b = [a0, a1];
262 GTnonlinOrder=DRNLnonlinearOrder;
263 GTnonlinBdry1=cell(nBFs, GTnonlinOrder);
264 GTnonlinBdry2=cell(nBFs, GTnonlinOrder);
265
266 % complete BM record (BM displacement)
267 DRNLoutput=zeros(nBFs, signalLength);
268
269
270 %% IHC ---
271 % IHC cilia activity and receptor potential
272 % viscous coupling between BM and stereocilia displacement
273 % Nyquist=sampleRate/2;
274 % IHCcutoff=1/(2*pi*IHC_cilia_RPParams.tc);
275 % [IHCciliaFilter_b,IHCciliaFilter_a]=...
276 % butter(1, IHCcutoff/Nyquist, 'high');
277 a1=dt/IHC_cilia_RPParams.tc-1; a0=1;
278 b0=1+ a1;
279 % high pass (i.e. low pass reversed)
280 IHCciliaFilter_b=[a0 a1]; IHCciliaFilter_a=b0;
281 % figure(9), freqz(IHCciliaFilter_b, IHCciliaFilter_a)
282
283 IHCciliaBndry=cell(nBFs,1);
284
285 % IHC apical conductance (Boltzman function)
286 IHC_C= IHC_cilia_RPParams.C;
287 IHCu0= IHC_cilia_RPParams.u0;
288 IHCu1= IHC_cilia_RPParams.u1;
289 IHCs0= IHC_cilia_RPParams.s0;
290 IHCs1= IHC_cilia_RPParams.s1;
291 IHCGmax= IHC_cilia_RPParams.Gmax;
292 IHCGa= IHC_cilia_RPParams.Ga; % (leakage)
293
294 IHCGu0 = IHCGa+IHCGmax./(1+exp(IHCu0/IHCs0).*(1+exp(IHCu1/IHCs1)));
295 IHCrestingCiliaCond=IHCGu0;
296
297 % Receptor potential
298 IHC_Cab= IHC_cilia_RPParams.Cab;
299 IHC_Gk= IHC_cilia_RPParams.Gk;
300 IHC_Et= IHC_cilia_RPParams.Et;
301 IHC_Ek= IHC_cilia_RPParams.Ek;
302 IHC_Ekp= IHC_Ek+IHC_Et*IHC_cilia_RPParams.Rpc;
303
304 IHCrestingV= (IHC_Gk*IHC_Ekp+IHCGu0*IHC_Et)/(IHCGu0+IHC_Gk);
305
306 IHC_Vnow= IHCrestingV*ones(nBFs,1); % initial voltage
307 IHC_RP= zeros(nBFs,segmentLength);
308
309 % complete record of IHC receptor potential (V)
310 IHCciliaDisplacement= zeros(nBFs,segmentLength);
311 IHCoutput= zeros(nBFs,signalLength);
312 IHC_cilia_output= zeros(nBFs,signalLength);
313
314
315 %% pre-synapse ---
316 % Each BF is replicated using a different fiber type to make a 'channel'
317 % The number of channels is nBFs x nANfiberTypes
318 % Fiber types are specified in terms of tauCa
319 nANfiberTypes= length(IHCpreSynapseParams.tauCa);
320 ANtauCas= IHCpreSynapseParams.tauCa;
321 nANchannels= nANfiberTypes*nBFs;
322 synapticCa= zeros(nANchannels,segmentLength);
323
324 % Calcium control (more calcium, greater release rate)
325 ECa=IHCpreSynapseParams.ECa;
326 gamma=IHCpreSynapseParams.gamma;
327 beta=IHCpreSynapseParams.beta;
328 tauM=IHCpreSynapseParams.tauM;
329 mICa=zeros(nANchannels,segmentLength);
330 GmaxCa=IHCpreSynapseParams.GmaxCa;
331 synapse_z= IHCpreSynapseParams.z;
332 synapse_power=IHCpreSynapseParams.power;
333
334 % tauCa vector is established across channels to allow vectorization
335 % (one tauCa per channel). Do not confuse with ANtauCas (one pre fiber type)
336 tauCa=repmat(ANtauCas, nBFs,1);
337 tauCa=reshape(tauCa, nANchannels, 1);
338
339 % presynapse startup values (vectors, length:nANchannels)
340 % proportion (0 - 1) of Ca channels open at IHCrestingV
341 mICaCurrent=((1+beta^-1 * exp(-gamma*IHCrestingV))^-1)...
342 *ones(nBFs*nANfiberTypes,1);
343 % corresponding startup currents
344 ICaCurrent= (GmaxCa*mICaCurrent.^3) * (IHCrestingV-ECa);
345 CaCurrent= ICaCurrent.*tauCa;
346
347 % vesicle release rate at startup (one per channel)
348 % kt0 is used only at initialisation
349 kt0= -synapse_z * CaCurrent.^synapse_power;
350
351
352 %% AN ---
353 % each row of the AN matrices represents one AN fiber
354 % The results computed either for probabiities *or* for spikes (not both)
355 % Spikes are necessary if CN and IC are to be computed
356 nFibersPerChannel= AN_IHCsynapseParams.numFibers;
357 nANfibers= nANchannels*nFibersPerChannel;
358 AN_refractory_period= AN_IHCsynapseParams.refractory_period;
359
360 y=AN_IHCsynapseParams.y;
361 l=AN_IHCsynapseParams.l;
362 x=AN_IHCsynapseParams.x;
363 r=AN_IHCsynapseParams.r;
364 M=round(AN_IHCsynapseParams.M);
365
366 % probability (NB initial 'P' on everything)
367 PAN_ydt = repmat(AN_IHCsynapseParams.y*dt, nANchannels,1);
368 PAN_ldt = repmat(AN_IHCsynapseParams.l*dt, nANchannels,1);
369 PAN_xdt = repmat(AN_IHCsynapseParams.x*dt, nANchannels,1);
370 PAN_rdt = repmat(AN_IHCsynapseParams.r*dt, nANchannels,1);
371 PAN_rdt_plus_ldt = PAN_rdt + PAN_ldt;
372 PAN_M=round(AN_IHCsynapseParams.M);
373
374 % compute starting values
375 Pcleft = kt0* y* M ./ (y*(l+r)+ kt0* l);
376 Pavailable = Pcleft*(l+r)./kt0;
377 Preprocess = Pcleft*r/x; % canbe fractional
378
379 ANprobability=zeros(nANchannels,segmentLength);
380 ANprobRateOutput=zeros(nANchannels,signalLength);
381 lengthAbsRefractoryP= round(AN_refractory_period/dt);
382 cumANnotFireProb=ones(nANchannels,signalLength);
383 % special variables for monitoring synaptic cleft (specialists only)
384 savePavailableSeg=zeros(nANchannels,segmentLength);
385 savePavailable=zeros(nANchannels,signalLength);
386
387 % spikes % ! ! ! ! ! ! ! !
388 lengthAbsRefractory= round(AN_refractory_period/ANdt);
389
390 AN_ydt= repmat(AN_IHCsynapseParams.y*ANdt, nANfibers,1);
391 AN_ldt= repmat(AN_IHCsynapseParams.l*ANdt, nANfibers,1);
392 AN_xdt= repmat(AN_IHCsynapseParams.x*ANdt, nANfibers,1);
393 AN_rdt= repmat(AN_IHCsynapseParams.r*ANdt, nANfibers,1);
394 AN_rdt_plus_ldt= AN_rdt + AN_ldt;
395 AN_M= round(AN_IHCsynapseParams.M);
396
397 % kt0 is initial release rate
398 % Establish as a vector (length=channel x number of fibers)
399 kt0= repmat(kt0', nFibersPerChannel, 1);
400 kt0=reshape(kt0, nANfibers,1);
401
402 % starting values for reservoirs
403 AN_cleft = kt0* y* M ./ (y*(l+r)+ kt0* l);
404 AN_available = round(AN_cleft*(l+r)./kt0); %must be integer
405 AN_reprocess = AN_cleft*r/x;
406
407 % output is in a logical array spikes = 1/ 0.
408 ANspikes= false(nANfibers,reducedSegmentLength);
409 ANoutput= false(nANfibers,reducedSignalLength);
410
411
412 %% CN (first brain stem nucleus - could be any subdivision of CN)
413 % Input to a CN neuorn is a random selection of AN fibers within a channel
414 % The number of AN fibers used is ANfibersFanInToCN
415 % CNtauGk (Potassium time constant) determines the rate of firing of
416 % the unit when driven hard by a DC input (not normally >350 sp/s)
417 % If there is more than one value, everything is replicated accordingly
418
419 ANavailableFibersPerChan=AN_IHCsynapseParams.numFibers;
420 ANfibersFanInToCN=MacGregorMultiParams.fibersPerNeuron;
421
422 CNtauGk=MacGregorMultiParams.tauGk; % row vector of CN types (by tauGk)
423 nCNtauGk=length(CNtauGk);
424
425 % the total number of 'channels' is now greater
426 % 'channel' is defined as collections of units with the same parameters
427 % i.e. same BF, same ANtau, same CNtauGk
428 nCNchannels=nANchannels*nCNtauGk;
429
430 nCNneuronsPerChannel=MacGregorMultiParams.nNeuronsPerBF;
431 tauGk=repmat(CNtauGk, nCNneuronsPerChannel,1);
432 tauGk=reshape(tauGk,nCNneuronsPerChannel*nCNtauGk,1);
433
434 % Now the number of neurons has been increased
435 nCNneurons=nCNneuronsPerChannel*nCNchannels;
436 CNmembranePotential=zeros(nCNneurons,reducedSegmentLength);
437
438 % establish which ANfibers (by name) feed into which CN nuerons
439 CNinputfiberLists=zeros(nANchannels*nCNneuronsPerChannel, ANfibersFanInToCN);
440 unitNo=1;
441 for ch=1:nANchannels
442 % Each channel contains a number of units =length(listOfFanInValues)
443 for idx=1:nCNneuronsPerChannel
444 for idx2=1:nCNtauGk
445 fibersUsed=(ch-1)*ANavailableFibersPerChan + ...
446 ceil(rand(1,ANfibersFanInToCN)* ANavailableFibersPerChan);
447 CNinputfiberLists(unitNo,:)=fibersUsed;
448 unitNo=unitNo+1;
449 end
450 end
451 end
452
453 % input to CN units
454 AN_PSTH=zeros(nCNneurons,reducedSegmentLength);
455
456 % Generate CNalphaFunction function
457 % by which spikes are converted to post-synaptic currents
458 CNdendriteLPfreq= MacGregorMultiParams.dendriteLPfreq;
459 CNcurrentPerSpike=MacGregorMultiParams.currentPerSpike;
460 CNspikeToCurrentTau=1/(2*pi*CNdendriteLPfreq);
461 t=ANdt:ANdt:5*CNspikeToCurrentTau;
462 CNalphaFunction= (1 / ...
463 CNspikeToCurrentTau)*t.*exp(-t /CNspikeToCurrentTau);
464 CNalphaFunction=CNalphaFunction*CNcurrentPerSpike;
465
466 % figure(98), plot(t,CNalphaFunction)
467 % working memory for implementing convolution
468
469 CNcurrentTemp=...
470 zeros(nCNneurons,reducedSegmentLength+length(CNalphaFunction)-1);
471 % trailing alphas are parts of humps carried forward to the next segment
472 CNtrailingAlphas=zeros(nCNneurons,length(CNalphaFunction));
473
474 CN_tauM=MacGregorMultiParams.tauM;
475 CN_tauTh=MacGregorMultiParams.tauTh;
476 CN_cap=MacGregorMultiParams.Cap;
477 CN_c=MacGregorMultiParams.c;
478 CN_b=MacGregorMultiParams.dGkSpike;
479 CN_Ek=MacGregorMultiParams.Ek;
480 CN_Eb= MacGregorMultiParams.Eb;
481 CN_Er=MacGregorMultiParams.Er;
482 CN_Th0= MacGregorMultiParams.Th0;
483 CN_E= zeros(nCNneurons,1);
484 CN_Gk= zeros(nCNneurons,1);
485 CN_Th= MacGregorMultiParams.Th0*ones(nCNneurons,1);
486 CN_Eb=CN_Eb.*ones(nCNneurons,1);
487 CN_Er=CN_Er.*ones(nCNneurons,1);
488 CNtimeSinceLastSpike=zeros(nCNneurons,1);
489 % tauGk is the main distinction between neurons
490 % in fact they are all the same in the standard model
491 tauGk=repmat(tauGk,nANchannels,1);
492
493 CNoutput=false(nCNneurons,reducedSignalLength);
494
495
496 %% MacGregor (IC - second nucleus) --------
497 nICcells=nANchannels*nCNtauGk; % one cell per channel
498 CN_PSTH=zeros(nICcells ,reducedSegmentLength);
499
500 % ICspikeWidth=0.00015; % this may need revisiting
501 % epochsPerSpike=round(ICspikeWidth/ANdt);
502 % if epochsPerSpike<1
503 % error(['MacGregorMulti: sample rate too low to support ' ...
504 % num2str(ICspikeWidth*1e6) ' microsec spikes']);
505 % end
506
507 % short names
508 IC_tauM=MacGregorParams.tauM;
509 IC_tauGk=MacGregorParams.tauGk;
510 IC_tauTh=MacGregorParams.tauTh;
511 IC_cap=MacGregorParams.Cap;
512 IC_c=MacGregorParams.c;
513 IC_b=MacGregorParams.dGkSpike;
514 IC_Th0=MacGregorParams.Th0;
515 IC_Ek=MacGregorParams.Ek;
516 IC_Eb= MacGregorParams.Eb;
517 IC_Er=MacGregorParams.Er;
518
519 IC_E=zeros(nICcells,1);
520 IC_Gk=zeros(nICcells,1);
521 IC_Th=IC_Th0*ones(nICcells,1);
522
523 % Dendritic filtering, all spikes are replaced by CNalphaFunction functions
524 ICdendriteLPfreq= MacGregorParams.dendriteLPfreq;
525 ICcurrentPerSpike=MacGregorParams.currentPerSpike;
526 ICspikeToCurrentTau=1/(2*pi*ICdendriteLPfreq);
527 t=ANdt:ANdt:3*ICspikeToCurrentTau;
528 IC_CNalphaFunction= (ICcurrentPerSpike / ...
529 ICspikeToCurrentTau)*t.*exp(-t / ICspikeToCurrentTau);
530 % figure(98), plot(t,IC_CNalphaFunction)
531
532 % working space for implementing alpha function
533 ICcurrentTemp=...
534 zeros(nICcells,reducedSegmentLength+length(IC_CNalphaFunction)-1);
535 ICtrailingAlphas=zeros(nICcells, length(IC_CNalphaFunction));
536
537 ICfiberTypeRates=zeros(nANfiberTypes,reducedSignalLength);
538 ICoutput=false(nICcells,reducedSignalLength);
539
540 ICmembranePotential=zeros(nICcells,reducedSegmentLength);
541 ICmembraneOutput=zeros(nICcells,signalLength);
542
543
544 %% Main program %% %% %% %% %% %% %% %% %% %% %% %% %% %%
545
546 % Compute the entire model for each segment
547 segmentStartPTR=1;
548 reducedSegmentPTR=1; % when sampling rate is reduced
549 while segmentStartPTR<signalLength
550 segmentEndPTR=segmentStartPTR+segmentLength-1;
551 % shorter segments after speed up.
552 shorterSegmentEndPTR=reducedSegmentPTR+reducedSegmentLength-1;
553
554 inputPressureSegment=inputSignal...
555 (:,segmentStartPTR:segmentStartPTR+segmentLength-1);
556
557 % segment debugging plots
558 % figure(98)
559 % plot(segmentTime,inputPressureSegment), title('signalSegment')
560
561
562 % OME ----------------------
563
564 % OME Stage 1: external resonances. Add to inputSignal pressure wave
565 y=inputPressureSegment;
566 for n=1:nOMEExtFilters
567 % any number of resonances can be used
568 [x OMEExtFilterBndry{n}] = ...
569 filter(ExtFilter_b{n},ExtFilter_a{n},...
570 inputPressureSegment, OMEExtFilterBndry{n});
571 x= x* OMEgainScalars(n);
572 % This is a parallel resonance so add it
573 y=y+x;
574 end
575 inputPressureSegment=y;
576 OMEextEarPressure(segmentStartPTR:segmentEndPTR)= inputPressureSegment;
577
578 % OME stage 2: convert input pressure (velocity) to
579 % tympanic membrane(TM) displacement using low pass filter
580 [TMdisplacementSegment OME_TMdisplacementBndry] = ...
581 filter(TMdisp_b,TMdisp_a,inputPressureSegment, ...
582 OME_TMdisplacementBndry);
583 % and save it
584 TMoutput(segmentStartPTR:segmentEndPTR)= TMdisplacementSegment;
585
586 % OME stage 3: middle ear high pass effect to simulate stapes inertia
587 [stapesDisplacement OMEhighPassBndry] = ...
588 filter(stapesDisp_b,stapesDisp_a,TMdisplacementSegment, ...
589 OMEhighPassBndry);
590
591 % OME stage 4: apply stapes scalar
592 stapesDisplacement=stapesDisplacement*stapesScalar;
593
594 % OME stage 5: acoustic reflex stapes attenuation
595 % Attenuate the TM response using feedback from LSR fiber activity
596 if segmentStartPTR>efferentDelayPts
597 stapesDisplacement= stapesDisplacement.*...
598 ARattenuation(segmentStartPTR-efferentDelayPts:...
599 segmentEndPTR-efferentDelayPts);
600 end
601
602 % segment debugging plots
603 % figure(98)
604 % plot(segmentTime, stapesDisplacement), title ('stapesDisplacement')
605
606 % and save
607 OMEoutput(segmentStartPTR:segmentEndPTR)= stapesDisplacement;
608
609
610 %% BM ------------------------------
611 % Each location is computed separately
612 for BFno=1:nBFs
613
614 % *linear* path
615 linOutput = stapesDisplacement * linGAIN; % linear gain
616
617 for order = 1 : GTlinOrder
618 [linOutput GTlinBdry{BFno,order}] = ...
619 filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, ...
620 GTlinBdry{BFno,order});
621 end
622
623 % *nonLinear* path
624 % efferent attenuation (0 <> 1)
625 if segmentStartPTR>efferentDelayPts
626 MOC=MOCattenuation(BFno, segmentStartPTR-efferentDelayPts:...
627 segmentEndPTR-efferentDelayPts);
628 else % no MOC available yet
629 MOC=ones(1, segmentLength);
630 end
631 % apply MOC to nonlinear input function
632 nonlinOutput=stapesDisplacement.* MOC;
633
634 % first gammatone filter (nonlin path)
635 for order = 1 : GTnonlinOrder
636 [nonlinOutput GTnonlinBdry1{BFno,order}] = ...
637 filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ...
638 nonlinOutput, GTnonlinBdry1{BFno,order});
639 end
640
641 % % original broken stick instantaneous compression
642 % y= nonlinOutput.* DRNLa; % linear section.
643 % % compress parts of the signal above the compression threshold
644 % abs_x = abs(nonlinOutput);
645 % idx=find(abs_x>DRNLcompressionThreshold);
646 % if ~isempty(idx)>0
647 % y(idx)=sign(y(idx)).* (DRNLb*abs_x(idx).^DRNLc);
648 % end
649 % nonlinOutput=y;
650
651
652 % new broken stick instantaneous compression
653 y= nonlinOutput.* DRNLa; % linear section attenuation/gain.
654 % compress parts of the signal above the compression threshold
655 % holdY=y;
656 abs_y = abs(y);
657 idx=find(abs_y>DRNLcompressionThreshold);
658 if ~isempty(idx)>0
659 % y(idx)=sign(y(idx)).* (DRNLcompressionThreshold + ...
660 % (abs_y(idx)-DRNLcompressionThreshold).^DRNLc);
661 y(idx)=sign(y(idx)).* (DRNLcompressionThreshold + ...
662 (abs_y(idx)-DRNLcompressionThreshold)*DRNLc);
663 end
664 nonlinOutput=y;
665
666 % % Boltzmann compression
667 % y=(nonlinOutput * DRNLa*100000);
668 % holdY=y;
669 % y=abs(y);
670 % s=10; u=0.0;
671 % x=1./(1+exp(-(y-u)/s))-0.5;
672 % nonlinOutput=sign(nonlinOutput).*x/10000;
673
674
675 % if segmentStartPTR==10*segmentLength+1
676 % figure(90)
677 % plot(holdY,'b'), hold on
678 % plot(nonlinOutput, 'r'), hold off
679 % ylim([-1e-5 1e-5])
680 % pause(1)
681 % end
682
683 % second filter removes distortion products
684 for order = 1 : GTnonlinOrder
685 [ nonlinOutput GTnonlinBdry2{BFno,order}] = ...
686 filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ...
687 nonlinOutput, GTnonlinBdry2{BFno,order});
688 end
689
690 % combine the two paths to give the DRNL displacement
691 DRNLresponse(BFno,:)=linOutput+nonlinOutput;
692 % disp(num2str(max(linOutput)))
693 end % BF
694
695 % segment debugging plots
696 % figure(98)
697 % if size(DRNLresponse,1)>3
698 % imagesc(DRNLresponse) % matrix display
699 % title('DRNLresponse'); % single or double channel response
700 % else
701 % plot(segmentTime, DRNLresponse)
702 % end
703
704 % and save it
705 DRNLoutput(:, segmentStartPTR:segmentEndPTR)= DRNLresponse;
706
707
708 %% IHC ------------------------------------
709 % BM displacement to IHCciliaDisplacement is a high-pass filter
710 % because of viscous coupling
711 for idx=1:nBFs
712 [IHCciliaDisplacement(idx,:) IHCciliaBndry{idx}] = ...
713 filter(IHCciliaFilter_b,IHCciliaFilter_a, ...
714 DRNLresponse(idx,:), IHCciliaBndry{idx});
715 end
716
717 % apply scalar
718 IHCciliaDisplacement=IHCciliaDisplacement* IHC_C;
719
720 % and save it
721 IHC_cilia_output(:,segmentStartPTR:segmentStartPTR+segmentLength-1)=...
722 IHCciliaDisplacement;
723
724 % compute apical conductance
725 G=IHCGmax./(1+exp(-(IHCciliaDisplacement-IHCu0)/IHCs0).*...
726 (1+exp(-(IHCciliaDisplacement-IHCu1)/IHCs1)));
727 Gu=G + IHCGa;
728
729 % Compute receptor potential
730 for idx=1:segmentLength
731 IHC_Vnow=IHC_Vnow+ (-Gu(:, idx).*(IHC_Vnow-IHC_Et)-...
732 IHC_Gk*(IHC_Vnow-IHC_Ekp))* dt/IHC_Cab;
733 IHC_RP(:,idx)=IHC_Vnow;
734 end
735
736 % segment debugging plots
737 % if size(IHC_RP,1)>3
738 % surf(IHC_RP), shading interp, title('IHC_RP')
739 % else
740 % plot(segmentTime, IHC_RP)
741 % end
742
743 % and save it
744 IHCoutput(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=IHC_RP;
745
746
747 %% synapse -----------------------------
748 % Compute the vesicle release rate for each fiber type at each BF
749 % replicate IHC_RP for each fiber type
750 Vsynapse=repmat(IHC_RP, nANfiberTypes,1);
751
752 % look-up table of target fraction channels open for a given IHC_RP
753 mICaINF= 1./( 1 + exp(-gamma * Vsynapse) /beta);
754 % fraction of channel open - apply time constant
755 for idx=1:segmentLength
756 % mICaINF is the current 'target' value of mICa
757 mICaCurrent=mICaCurrent+(mICaINF(:,idx)-mICaCurrent)*dt./tauM;
758 mICa(:,idx)=mICaCurrent;
759 end
760
761 ICa= (GmaxCa* mICa.^3) .* (Vsynapse- ECa);
762
763 for idx=1:segmentLength
764 CaCurrent=CaCurrent + ICa(:,idx)*dt - CaCurrent*dt./tauCa;
765 synapticCa(:,idx)=CaCurrent;
766 end
767 synapticCa=-synapticCa; % treat synapticCa as positive substance
768
769 % NB vesicleReleaseRate is /s and is independent of dt
770 vesicleReleaseRate = synapse_z * synapticCa.^synapse_power; % rate
771
772 % segment debugging plots
773 % if size(vesicleReleaseRate,1)>3
774 % surf(vesicleReleaseRate), shading interp, title('vesicleReleaseRate')
775 % else
776 % plot(segmentTime, vesicleReleaseRate)
777 % end
778
779
780 %% AN
781 switch AN_spikesOrProbability
782 case 'probability'
783 % No refractory effect is applied
784 for t = 1:segmentLength;
785 M_Pq=PAN_M-Pavailable;
786 M_Pq(M_Pq<0)=0;
787 Preplenish = M_Pq .* PAN_ydt;
788 Pejected = Pavailable.* vesicleReleaseRate(:,t)*dt;
789 Preprocessed = M_Pq.*Preprocess.* PAN_xdt;
790
791 ANprobability(:,t)= min(Pejected,1);
792 reuptakeandlost= PAN_rdt_plus_ldt .* Pcleft;
793 reuptake= PAN_rdt.* Pcleft;
794
795 Pavailable= Pavailable+ Preplenish- Pejected+ Preprocessed;
796 Pcleft= Pcleft + Pejected - reuptakeandlost;
797 Preprocess= Preprocess + reuptake - Preprocessed;
798 Pavailable(Pavailable<0)=0;
799 savePavailableSeg(:,t)=Pavailable; % synapse tracking
800
801 end
802
803 % and save it as *rate*
804 ANrate=ANprobability/dt;
805 ANprobRateOutput(:, segmentStartPTR:...
806 segmentStartPTR+segmentLength-1)= ANrate;
807 % monitor synapse contents (only sometimes used)
808 savePavailable(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=...
809 savePavailableSeg;
810
811 %% Apply refractory effect
812 % the probability of a spike's occurring in the preceding
813 % refractory window (t= tnow-refractory period to tnow-)
814 % pFired= 1 - II(1-p(t)),
815 % we need a running account of cumProb=II(1-p(t))
816 % cumProb(t)= cumProb(t-1)*(1-p(t))/(1-p(t-refracPeriod))
817 % cumProb(0)=0
818 % pFired(t)= 1-cumProb(t)
819 % This gives the fraction of firing events that must be
820 % discounted because of a firing event in the refractory
821 % period
822 % p(t)= ANprobOutput(t) * pFired(t)
823 % where ANprobOutput is the uncorrected firing probability
824 % based on vesicle release rate
825 % NB this covers only the absoute refractory period
826 % not the relative refractory period. To approximate this it
827 % is necessary to extend the refractory period by 50%
828
829
830 for t = segmentStartPTR:segmentEndPTR;
831 if t>1
832 ANprobRateOutput(:,t)= ANprobRateOutput(:,t)...
833 .* cumANnotFireProb(:,t-1);
834 end
835 % add recent and remove distant probabilities
836 refrac=round(lengthAbsRefractoryP * 1.5);
837 if t>refrac
838 cumANnotFireProb(:,t)= cumANnotFireProb(:,t-1)...
839 .*(1-ANprobRateOutput(:,t)*dt)...
840 ./(1-ANprobRateOutput(:,t-refrac)*dt);
841 end
842 end
843 % figure(88), plot(cumANnotFireProb'), title('cumNotFire')
844 % figure(89), plot(ANprobRateOutput'), title('ANprobRateOutput')
845
846 %% Estimate efferent effects. ARattenuation (0 <> 1)
847 % acoustic reflex
848 [r c]=size(ANrate);
849 if r>nBFs % Only if LSR fibers are computed
850 ARAttSeg=mean(ANrate(1:nBFs,:),1); %LSR channels are 1:nBF
851 % smooth
852 [ARAttSeg, ARboundaryProb] = ...
853 filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundaryProb);
854 ARAttSeg=ARAttSeg-ARrateThreshold;
855 ARAttSeg(ARAttSeg<0)=0; % prevent negative strengths
856 ARattenuation(segmentStartPTR:segmentEndPTR)=...
857 (1-ARrateToAttenuationFactorProb.* ARAttSeg);
858 end
859 % plot(ARattenuation)
860
861 % MOC attenuation based on within-channel HSR fiber activity
862 HSRbegins=nBFs*(nANfiberTypes-1)+1;
863 rates=ANrate(HSRbegins:end,:);
864 if rateToAttenuationFactorProb<0
865 % negative factor implies a fixed attenuation
866 MOCattenuation(:,segmentStartPTR:segmentEndPTR)= ...
867 ones(size(rates))* -rateToAttenuationFactorProb;
868 else
869 for idx=1:nBFs
870 [smoothedRates, MOCprobBoundary{idx}] = ...
871 filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ...
872 MOCprobBoundary{idx});
873 smoothedRates=smoothedRates-MOCrateThresholdProb;
874 smoothedRates=max(smoothedRates, 0);
875
876 x=(1- smoothedRates* rateToAttenuationFactorProb);
877 x=max(x, 10^(-30/20));
878 MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= x;
879 end
880 end
881 MOCattenuation(MOCattenuation<0)=0.001;
882
883 % plot(MOCattenuation)
884
885
886 case 'spikes'
887 ANtimeCount=0;
888 % implement speed upt
889 for t = ANspeedUpFactor:ANspeedUpFactor:segmentLength;
890 ANtimeCount=ANtimeCount+1;
891 % convert release rate to probabilities
892 releaseProb=vesicleReleaseRate(:,t)*ANdt;
893 % releaseProb is the release probability per channel
894 % but each channel has many synapses
895 releaseProb=repmat(releaseProb',nFibersPerChannel,1);
896 releaseProb=reshape(releaseProb, nFibersPerChannel*nANchannels,1);
897
898 % AN_available=round(AN_available); % vesicles must be integer, (?needed)
899 M_q=AN_M- AN_available; % number of missing vesicles
900 M_q(M_q<0)= 0; % cannot be less than 0
901
902 % AN_N1 converts probability to discrete events
903 % it considers each event that might occur
904 % (how many vesicles might be released)
905 % and returns a count of how many were released
906
907 % slow line
908 % probabilities= 1-(1-releaseProb).^AN_available;
909 probabilities= 1-intpow((1-releaseProb), AN_available);
910 ejected= probabilities> rand(length(AN_available),1);
911
912 reuptakeandlost = AN_rdt_plus_ldt .* AN_cleft;
913 reuptake = AN_rdt.* AN_cleft;
914
915 % slow line
916 % probabilities= 1-(1-AN_reprocess.*AN_xdt).^M_q;
917 probabilities= 1-intpow((1-AN_reprocess.*AN_xdt), M_q);
918 reprocessed= probabilities>rand(length(M_q),1);
919
920 % slow line
921 % probabilities= 1-(1-AN_ydt).^M_q;
922 probabilities= 1-intpow((1-AN_ydt), M_q);
923
924 replenish= probabilities>rand(length(M_q),1);
925
926 AN_available = AN_available + replenish - ejected ...
927 + reprocessed;
928 AN_cleft = AN_cleft + ejected - reuptakeandlost;
929 AN_reprocess = AN_reprocess + reuptake - reprocessed;
930
931 % ANspikes is logical record of vesicle release events>0
932 ANspikes(:, ANtimeCount)= ejected;
933 end % t
934
935 % zero any events that are preceded by release events ...
936 % within the refractory period
937 % The refractory period consist of two periods
938 % 1) the absolute period where no spikes occur
939 % 2) a relative period where a spike may occur. This relative
940 % period is realised as a variable length interval
941 % where the length is chosen at random
942 % (esentially a linear ramp up)
943
944 % Andreas has a fix for this
945 for t = 1:ANtimeCount-2*lengthAbsRefractory;
946 % identify all spikes across fiber array at time (t)
947 % idx is a list of channels where spikes occurred
948 % ?? try sparse matrices?
949 idx=find(ANspikes(:,t));
950 for j=idx % consider each spike
951 % specify variable refractory period
952 % between abs and 2*abs refractory period
953 nPointsRefractory=lengthAbsRefractory+...
954 round(rand*lengthAbsRefractory);
955 % disable spike potential for refractory period
956 % set all values in this range to 0
957 ANspikes(j,t+1:t+nPointsRefractory)=0;
958 end
959 end %t
960
961 % segment debugging
962 % plotInstructions.figureNo=98;
963 % plotInstructions.displaydt=ANdt;
964 % plotInstructions.numPlots=1;
965 % plotInstructions.subPlotNo=1;
966 % UTIL_plotMatrix(ANspikes, plotInstructions);
967
968 % and save it. NB, AN is now on 'speedUp' time
969 ANoutput(:, reducedSegmentPTR: shorterSegmentEndPTR)=ANspikes;
970
971
972 %% CN Macgregor first neucleus -------------------------------
973 % input is from AN so ANdt is used throughout
974 % Each CNneuron has a unique set of input fibers selected
975 % at random from the available AN fibers (CNinputfiberLists)
976
977 % Create the dendritic current for that neuron
978 % First get input spikes to this neuron
979 synapseNo=1;
980 for ch=1:nCNchannels
981 for idx=1:nCNneuronsPerChannel
982 % determine candidate fibers for this unit
983 fibersUsed=CNinputfiberLists(synapseNo,:);
984 % ANpsth has a bin width of ANdt
985 % (just a simple sum across fibers)
986 AN_PSTH(synapseNo,:) = ...
987 sum(ANspikes(fibersUsed,:), 1);
988 synapseNo=synapseNo+1;
989 end
990 end
991
992 % One alpha function per spike
993 [alphaRows alphaCols]=size(CNtrailingAlphas);
994
995 for unitNo=1:nCNneurons
996 CNcurrentTemp(unitNo,:)= ...
997 conv2(AN_PSTH(unitNo,:),CNalphaFunction);
998 end
999 % disp(['sum(AN_PSTH)= ' num2str(sum(AN_PSTH(1,:)))])
1000 % add post-synaptic current left over from previous segment
1001 CNcurrentTemp(:,1:alphaCols)=...
1002 CNcurrentTemp(:,1:alphaCols)+ CNtrailingAlphas;
1003
1004 % take post-synaptic current for this segment
1005 CNcurrentInput= CNcurrentTemp(:, 1:reducedSegmentLength);
1006 % disp(['mean(CNcurrentInput)= ' num2str(mean(CNcurrentInput(1,:)))])
1007
1008 % trailingalphas are the ends of the alpha functions that
1009 % spill over into the next segment
1010 CNtrailingAlphas= ...
1011 CNcurrentTemp(:, reducedSegmentLength+1:end);
1012
1013 if CN_c>0
1014 % variable threshold condition (slow)
1015 for t=1:reducedSegmentLength
1016 CNtimeSinceLastSpike=CNtimeSinceLastSpike-ANdt;
1017 s=CN_E>CN_Th & CNtimeSinceLastSpike<0 ;
1018 CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike
1019 dE =(-CN_E/CN_tauM + ...
1020 CNcurrentInput(:,t)/CN_cap+(...
1021 CN_Gk/CN_cap).*(CN_Ek-CN_E))*ANdt;
1022 dGk=-CN_Gk*ANdt./tauGk + CN_b*s;
1023 dTh=-(CN_Th-CN_Th0)*ANdt/CN_tauTh + CN_c*s;
1024 CN_E=CN_E+dE;
1025 CN_Gk=CN_Gk+dGk;
1026 CN_Th=CN_Th+dTh;
1027 CNmembranePotential(:,t)=CN_E+s.*(CN_Eb-CN_E)+CN_Er;
1028 end
1029 else
1030 % static threshold (faster)
1031 E=zeros(1,reducedSegmentLength);
1032 Gk=zeros(1,reducedSegmentLength);
1033 ss=zeros(1,reducedSegmentLength);
1034 for t=1:reducedSegmentLength
1035 % time of previous spike moves back in time
1036 CNtimeSinceLastSpike=CNtimeSinceLastSpike-ANdt;
1037 % action potential if E>threshold
1038 % allow time for s to reset between events
1039 s=CN_E>CN_Th0 & CNtimeSinceLastSpike<0 ;
1040 ss(t)=s(1);
1041 CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike
1042 dE = (-CN_E/CN_tauM + ...
1043 CNcurrentInput(:,t)/CN_cap +...
1044 (CN_Gk/CN_cap).*(CN_Ek-CN_E))*ANdt;
1045 dGk=-CN_Gk*ANdt./tauGk +CN_b*s;
1046 CN_E=CN_E+dE;
1047 CN_Gk=CN_Gk+dGk;
1048 E(t)=CN_E(1);
1049 Gk(t)=CN_Gk(1);
1050 % add spike to CN_E and add resting potential (-60 mV)
1051 CNmembranePotential(:,t)=CN_E +s.*(CN_Eb-CN_E)+CN_Er;
1052 end
1053 end
1054 % disp(['CN_E= ' num2str(sum(CN_E(1,:)))])
1055 % disp(['CN_Gk= ' num2str(sum(CN_Gk(1,:)))])
1056 % disp(['CNmembranePotential= ' num2str(sum(CNmembranePotential(1,:)))])
1057 % plot(CNmembranePotential(1,:))
1058
1059
1060 % extract spikes. A spike is a substantial upswing in voltage
1061 CN_spikes=CNmembranePotential> -0.02;
1062 % disp(['CNspikesbefore= ' num2str(sum(sum(CN_spikes)))])
1063
1064 % now remove any spike that is immediately followed by a spike
1065 % NB 'find' works on columns (whence the transposing)
1066 % for each spike put a zero in the next epoch
1067 CN_spikes=CN_spikes';
1068 idx=find(CN_spikes);
1069 idx=idx(1:end-1);
1070 CN_spikes(idx+1)=0;
1071 CN_spikes=CN_spikes';
1072 % disp(['CNspikes= ' num2str(sum(sum(CN_spikes)))])
1073
1074 % segment debugging
1075 % plotInstructions.figureNo=98;
1076 % plotInstructions.displaydt=ANdt;
1077 % plotInstructions.numPlots=1;
1078 % plotInstructions.subPlotNo=1;
1079 % UTIL_plotMatrix(CN_spikes, plotInstructions);
1080
1081 % and save it
1082 CNoutput(:, reducedSegmentPTR:shorterSegmentEndPTR)=...
1083 CN_spikes;
1084
1085
1086 %% IC ----------------------------------------------
1087 % MacGregor or some other second order neurons
1088
1089 % combine CN neurons in same channel (BF x AN tau x CNtau)
1090 % i.e. same BF, same tauCa, same CNtau
1091 % to generate inputs to single IC unit
1092 channelNo=0;
1093 for idx=1:nCNneuronsPerChannel: ...
1094 nCNneurons-nCNneuronsPerChannel+1;
1095 channelNo=channelNo+1;
1096 CN_PSTH(channelNo,:)=...
1097 sum(CN_spikes(idx:idx+nCNneuronsPerChannel-1,:));
1098 end
1099
1100 [alphaRows alphaCols]=size(ICtrailingAlphas);
1101 for ICneuronNo=1:nICcells
1102 ICcurrentTemp(ICneuronNo,:)= ...
1103 conv2(CN_PSTH(ICneuronNo,:), IC_CNalphaFunction);
1104 end
1105
1106 % add the unused current from the previous convolution
1107 ICcurrentTemp(:,1:alphaCols)=ICcurrentTemp(:,1:alphaCols)...
1108 + ICtrailingAlphas;
1109 % take what is required and keep the trailing part for next time
1110 inputCurrent=ICcurrentTemp(:, 1:reducedSegmentLength);
1111 ICtrailingAlphas=ICcurrentTemp(:, reducedSegmentLength+1:end);
1112
1113 if IC_c==0
1114 % faster computation when threshold is stable (c==0)
1115 for t=1:reducedSegmentLength
1116 s=IC_E>IC_Th0;
1117 dE = (-IC_E/IC_tauM + inputCurrent(:,t)/IC_cap +...
1118 (IC_Gk/IC_cap).*(IC_Ek-IC_E))*ANdt;
1119 dGk=-IC_Gk*ANdt/IC_tauGk +IC_b*s;
1120 IC_E=IC_E+dE;
1121 IC_Gk=IC_Gk+dGk;
1122 ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er;
1123 end
1124 else
1125 % threshold is changing (IC_c>0; e.g. bushy cell)
1126 for t=1:reducedSegmentLength
1127 dE = (-IC_E/IC_tauM + ...
1128 inputCurrent(:,t)/IC_cap + (IC_Gk/IC_cap)...
1129 .*(IC_Ek-IC_E))*ANdt;
1130 IC_E=IC_E+dE;
1131 s=IC_E>IC_Th;
1132 ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er;
1133 dGk=-IC_Gk*ANdt/IC_tauGk +IC_b*s;
1134 IC_Gk=IC_Gk+dGk;
1135
1136 % After a spike, the threshold is raised
1137 % otherwise it settles to its baseline
1138 dTh=-(IC_Th-Th0)*ANdt/IC_tauTh +s*IC_c;
1139 IC_Th=IC_Th+dTh;
1140 end
1141 end
1142
1143 ICspikes=ICmembranePotential> -0.01;
1144 % now remove any spike that is immediately followed by a spike
1145 % NB 'find' works on columns (whence the transposing)
1146 ICspikes=ICspikes';
1147 idx=find(ICspikes);
1148 idx=idx(1:end-1);
1149 ICspikes(idx+1)=0;
1150 ICspikes=ICspikes';
1151
1152 nCellsPerTau= nICcells/nANfiberTypes;
1153 firstCell=1;
1154 lastCell=nCellsPerTau;
1155 for tauCount=1:nANfiberTypes
1156 % separate rates according to fiber types
1157 % currently only the last segment is saved
1158 ICfiberTypeRates(tauCount, ...
1159 reducedSegmentPTR:shorterSegmentEndPTR)=...
1160 sum(ICspikes(firstCell:lastCell, :))...
1161 /(nCellsPerTau*ANdt);
1162 firstCell=firstCell+nCellsPerTau;
1163 lastCell=lastCell+nCellsPerTau;
1164 end
1165
1166 ICoutput(:,reducedSegmentPTR:shorterSegmentEndPTR)=ICspikes;
1167
1168 % store membrane output on original dt scale
1169 % do this for single channel models only
1170 % and only for the HSR-driven IC cells
1171 if round(nICcells/length(ANtauCas))==1 % single channel
1172 % select HSR driven cells
1173 x= ICmembranePotential(length(ANtauCas),:);
1174 % restore original dt
1175 x= repmat(x, ANspeedUpFactor,1);
1176 x= reshape(x,1,segmentLength);
1177 if nANfiberTypes>1 % save HSR and LSR
1178 y=repmat(ICmembranePotential(end,:),...
1179 ANspeedUpFactor,1);
1180 y= reshape(y,1,segmentLength);
1181 x=[x; y];
1182 end
1183 ICmembraneOutput(:, segmentStartPTR:segmentEndPTR)= x;
1184 end
1185
1186 % estimate efferent effects.
1187 % ARis based on LSR units. LSR channels are 1:nBF
1188 if nANfiberTypes>1 % AR is multi-channel only
1189 ARAttSeg=sum(ICspikes(1:nBFs,:),1)/ANdt;
1190 [ARAttSeg, ARboundary] = ...
1191 filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundary);
1192 ARAttSeg=ARAttSeg-ARrateThreshold;
1193 ARAttSeg(ARAttSeg<0)=0; % prevent negative strengths
1194 % scale up to dt from ANdt
1195 x= repmat(ARAttSeg, ANspeedUpFactor,1);
1196 x=reshape(x,1,segmentLength);
1197 ARattenuation(segmentStartPTR:segmentEndPTR)=...
1198 (1-ARrateToAttenuationFactor* x);
1199 ARattenuation(ARattenuation<0)=0.001;
1200 else
1201 % single channel model; disable AR
1202 ARattenuation(segmentStartPTR:segmentEndPTR)=...
1203 ones(1,segmentLength);
1204 end
1205
1206 % MOC attenuation using HSR response only
1207 % Separate MOC effect for each BF
1208 HSRbegins=nBFs*(nANfiberTypes-1)+1;
1209 rates=ICspikes(HSRbegins:end,:)/ANdt;
1210 for idx=1:nBFs
1211 [smoothedRates, MOCboundary{idx}] = ...
1212 filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ...
1213 MOCboundary{idx});
1214 % spont 'rates' is zero for IC
1215 MOCattSegment(idx,:)=smoothedRates;
1216 % expand timescale back to model dt from ANdt
1217 x= repmat(MOCattSegment(idx,:), ANspeedUpFactor,1);
1218 x= reshape(x,1,segmentLength);
1219 MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ...
1220 (1- MOCrateToAttenuationFactor* x);
1221 end
1222 MOCattenuation(MOCattenuation<0)=0.04;
1223 % segment debugging
1224 % plotInstructions.figureNo=98;
1225 % plotInstructions.displaydt=ANdt;
1226 % plotInstructions.numPlots=1;
1227 % plotInstructions.subPlotNo=1;
1228 % UTIL_plotMatrix(ICspikes, plotInstructions);
1229
1230 end % AN_spikesOrProbability
1231 segmentStartPTR=segmentStartPTR+segmentLength;
1232 reducedSegmentPTR=reducedSegmentPTR+reducedSegmentLength;
1233
1234
1235 end % segment
1236
1237 %% apply refractory correction to spike probabilities
1238
1239 % the probability of a spike's having occurred in the preceding
1240 % refractory window
1241 % pFired= 1 - II(1-p(t)), t= tnow-refractory period: tnow-1
1242 % we need a running account of cumProb=II(1-p(t))
1243 % cumProb(t)= cumProb(t-1)*(1-p(t-1))/(1-p(t-refracPeriod))
1244 % cumProb(0)=0
1245 % pFired(t)= 1-cumProb(t)
1246 % whence
1247 % p(t)= ANprobOutput(t) * pFired(t)
1248 % where ANprobOutput is the uncorrected probability
1249
1250
1251 % switch AN_spikesOrProbability
1252 % case 'probability'
1253 % ANprobOutput=ANprobRateOutput*dt;
1254 % [r nEpochs]=size(ANprobOutput);
1255 % % find probability of no spikes in refractory period
1256 % pNoSpikesInRefrac=ones(size(ANprobOutput));
1257 % pSpike=zeros(size(ANprobOutput));
1258 % for epochNo=lengthAbsRefractoryP+2:nEpochs
1259 % pNoSpikesInRefrac(:,epochNo)=...
1260 % pNoSpikesInRefrac(:,epochNo-2)...
1261 % .*(1-pSpike(:,epochNo-1))...
1262 % ./(1-pSpike(:,epochNo-lengthAbsRefractoryP-1));
1263 % pSpike(:,epochNo)= ANprobOutput(:,epochNo)...
1264 % .*pNoSpikesInRefrac(:,epochNo);
1265 % end
1266 % ANprobRateOutput=pSpike/dt;
1267 % end
1268
1269 path(restorePath)