Revision 18:e9e263e4fcde MAP/MAP1_14.m
| MAP/MAP1_14.m | ||
|---|---|---|
| 83 | 83 |
segmentTime=dt*(1:segmentLength); % used in debugging plots |
| 84 | 84 |
|
| 85 | 85 |
% all spiking activity is computed using longer epochs |
| 86 |
ANspeedUpFactor=5; % 5 times longer
|
|
| 86 |
ANspeedUpFactor=AN_IHCsynapseParams.ANspeedUpFactor; % e.g.5 times
|
|
| 87 | 87 |
|
| 88 | 88 |
% inputSignal must be row vector |
| 89 | 89 |
[r c]=size(inputSignal); |
| ... | ... | |
| 443 | 443 |
CNcurrentPerSpike=MacGregorMultiParams.currentPerSpike; |
| 444 | 444 |
CNspikeToCurrentTau=1/(2*pi*CNdendriteLPfreq); |
| 445 | 445 |
t=ANdt:ANdt:5*CNspikeToCurrentTau; |
| 446 |
CNalphaFunction=... |
|
| 447 |
(CNcurrentPerSpike/CNspikeToCurrentTau)*t.*exp(-t/CNspikeToCurrentTau); |
|
| 446 |
CNalphaFunction= (1 / ... |
|
| 447 |
CNspikeToCurrentTau)*t.*exp(-t /CNspikeToCurrentTau); |
|
| 448 |
CNalphaFunction=CNalphaFunction*CNcurrentPerSpike; |
|
| 449 |
|
|
| 448 | 450 |
% figure(98), plot(t,CNalphaFunction) |
| 449 | 451 |
% working memory for implementing convolution |
| 452 |
|
|
| 450 | 453 |
CNcurrentTemp=... |
| 451 | 454 |
zeros(nCNneurons,reducedSegmentLength+length(CNalphaFunction)-1); |
| 452 | 455 |
% trailing alphas are parts of humps carried forward to the next segment |
| ... | ... | |
| 879 | 882 |
for idx=1:nCNneuronsPerChannel |
| 880 | 883 |
% determine candidate fibers for this unit |
| 881 | 884 |
fibersUsed=CNinputfiberLists(synapseNo,:); |
| 882 |
% ANpsth has a bin width of dt |
|
| 885 |
% ANpsth has a bin width of ANdt
|
|
| 883 | 886 |
% (just a simple sum across fibers) |
| 884 | 887 |
AN_PSTH(synapseNo,:) = ... |
| 885 | 888 |
sum(ANspikes(fibersUsed,:), 1); |
| ... | ... | |
| 894 | 897 |
CNcurrentTemp(unitNo,:)= ... |
| 895 | 898 |
conv(AN_PSTH(unitNo,:),CNalphaFunction); |
| 896 | 899 |
end |
| 900 |
% disp(['sum(AN_PSTH)= ' num2str(sum(AN_PSTH(1,:)))]) |
|
| 897 | 901 |
% add post-synaptic current left over from previous segment |
| 898 | 902 |
CNcurrentTemp(:,1:alphaCols)=... |
| 899 | 903 |
CNcurrentTemp(:,1:alphaCols)+ CNtrailingAlphas; |
| 900 | 904 |
|
| 901 | 905 |
% take post-synaptic current for this segment |
| 902 | 906 |
CNcurrentInput= CNcurrentTemp(:, 1:reducedSegmentLength); |
| 907 |
% disp(['mean(CNcurrentInput)= ' num2str(mean(CNcurrentInput(1,:)))]) |
|
| 903 | 908 |
|
| 904 | 909 |
% trailingalphas are the ends of the alpha functions that |
| 905 | 910 |
% spill over into the next segment |
| ... | ... | |
| 909 | 914 |
if CN_c>0 |
| 910 | 915 |
% variable threshold condition (slow) |
| 911 | 916 |
for t=1:reducedSegmentLength |
| 912 |
CNtimeSinceLastSpike=CNtimeSinceLastSpike-dts;
|
|
| 917 |
CNtimeSinceLastSpike=CNtimeSinceLastSpike-ANdt;
|
|
| 913 | 918 |
s=CN_E>CN_Th & CNtimeSinceLastSpike<0 ; |
| 914 | 919 |
CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike |
| 915 | 920 |
dE =(-CN_E/CN_tauM + ... |
| 916 |
CNcurrentInput(:,t)/CN_cap+(CN_Gk/CN_cap).*(CN_Ek-CN_E))*dt; |
|
| 917 |
dGk=-CN_Gk*dt./tauGk + CN_b*s; |
|
| 918 |
dTh=-(CN_Th-CN_Th0)*dt/CN_tauTh + CN_c*s; |
|
| 921 |
CNcurrentInput(:,t)/CN_cap+(... |
|
| 922 |
CN_Gk/CN_cap).*(CN_Ek-CN_E))*ANdt; |
|
| 923 |
dGk=-CN_Gk*ANdt./tauGk + CN_b*s; |
|
| 924 |
dTh=-(CN_Th-CN_Th0)*ANdt/CN_tauTh + CN_c*s; |
|
| 919 | 925 |
CN_E=CN_E+dE; |
| 920 | 926 |
CN_Gk=CN_Gk+dGk; |
| 921 | 927 |
CN_Th=CN_Th+dTh; |
| ... | ... | |
| 923 | 929 |
end |
| 924 | 930 |
else |
| 925 | 931 |
% static threshold (faster) |
| 932 |
E=zeros(1,reducedSegmentLength); |
|
| 933 |
Gk=zeros(1,reducedSegmentLength); |
|
| 934 |
ss=zeros(1,reducedSegmentLength); |
|
| 926 | 935 |
for t=1:reducedSegmentLength |
| 927 |
CNtimeSinceLastSpike=CNtimeSinceLastSpike-dt; |
|
| 928 |
s=CN_E>CN_Th0 & CNtimeSinceLastSpike<0 ; % =1 if both conditions met |
|
| 929 |
CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike |
|
| 936 |
% time of previous spike moves back in time |
|
| 937 |
CNtimeSinceLastSpike=CNtimeSinceLastSpike-ANdt; |
|
| 938 |
% action potential if E>threshold |
|
| 939 |
% allow time for s to reset between events |
|
| 940 |
s=CN_E>CN_Th0 & CNtimeSinceLastSpike<0 ; |
|
| 941 |
ss(t)=s(1); |
|
| 942 |
CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike |
|
| 930 | 943 |
dE = (-CN_E/CN_tauM + ... |
| 931 |
CNcurrentInput(:,t)/CN_cap+(CN_Gk/CN_cap).*(CN_Ek-CN_E))*dt; |
|
| 932 |
dGk=-CN_Gk*dt./tauGk +CN_b*s; |
|
| 944 |
CNcurrentInput(:,t)/CN_cap +... |
|
| 945 |
(CN_Gk/CN_cap).*(CN_Ek-CN_E))*ANdt; |
|
| 946 |
dGk=-CN_Gk*ANdt./tauGk +CN_b*s; |
|
| 933 | 947 |
CN_E=CN_E+dE; |
| 934 | 948 |
CN_Gk=CN_Gk+dGk; |
| 949 |
E(t)=CN_E(1); |
|
| 950 |
Gk(t)=CN_Gk(1); |
|
| 935 | 951 |
% add spike to CN_E and add resting potential (-60 mV) |
| 936 |
CNmembranePotential(:,t)=CN_E+s.*(CN_Eb-CN_E)+CN_Er; |
|
| 952 |
CNmembranePotential(:,t)=CN_E +s.*(CN_Eb-CN_E)+CN_Er;
|
|
| 937 | 953 |
end |
| 938 | 954 |
end |
| 955 |
% disp(['CN_E= ' num2str(sum(CN_E(1,:)))]) |
|
| 956 |
% disp(['CN_Gk= ' num2str(sum(CN_Gk(1,:)))]) |
|
| 957 |
% disp(['CNmembranePotential= ' num2str(sum(CNmembranePotential(1,:)))]) |
|
| 958 |
% plot(CNmembranePotential(1,:)) |
|
| 959 |
|
|
| 939 | 960 |
|
| 940 | 961 |
% extract spikes. A spike is a substantial upswing in voltage |
| 941 |
CN_spikes=CNmembranePotential> -0.01; |
|
| 962 |
CN_spikes=CNmembranePotential> -0.02; |
|
| 963 |
% disp(['CNspikesbefore= ' num2str(sum(sum(CN_spikes)))]) |
|
| 942 | 964 |
|
| 943 | 965 |
% now remove any spike that is immediately followed by a spike |
| 944 | 966 |
% NB 'find' works on columns (whence the transposing) |
| 967 |
% for each spike put a zero in the next epoch |
|
| 945 | 968 |
CN_spikes=CN_spikes'; |
| 946 | 969 |
idx=find(CN_spikes); |
| 947 | 970 |
idx=idx(1:end-1); |
| 948 | 971 |
CN_spikes(idx+1)=0; |
| 949 | 972 |
CN_spikes=CN_spikes'; |
| 973 |
% disp(['CNspikes= ' num2str(sum(sum(CN_spikes)))]) |
|
| 950 | 974 |
|
| 951 | 975 |
% segment debugging |
| 952 | 976 |
% plotInstructions.figureNo=98; |
| ... | ... | |
| 990 | 1014 |
for t=1:reducedSegmentLength |
| 991 | 1015 |
s=IC_E>IC_Th0; |
| 992 | 1016 |
dE = (-IC_E/IC_tauM + inputCurrent(:,t)/IC_cap +... |
| 993 |
(IC_Gk/IC_cap).*(IC_Ek-IC_E))*dt; |
|
| 994 |
dGk=-IC_Gk*dt/IC_tauGk +IC_b*s; |
|
| 1017 |
(IC_Gk/IC_cap).*(IC_Ek-IC_E))*ANdt;
|
|
| 1018 |
dGk=-IC_Gk*ANdt/IC_tauGk +IC_b*s;
|
|
| 995 | 1019 |
IC_E=IC_E+dE; |
| 996 | 1020 |
IC_Gk=IC_Gk+dGk; |
| 997 | 1021 |
ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er; |
| ... | ... | |
| 1001 | 1025 |
for t=1:reducedSegmentLength |
| 1002 | 1026 |
dE = (-IC_E/IC_tauM + ... |
| 1003 | 1027 |
inputCurrent(:,t)/IC_cap + (IC_Gk/IC_cap)... |
| 1004 |
.*(IC_Ek-IC_E))*dt; |
|
| 1028 |
.*(IC_Ek-IC_E))*ANdt;
|
|
| 1005 | 1029 |
IC_E=IC_E+dE; |
| 1006 | 1030 |
s=IC_E>IC_Th; |
| 1007 | 1031 |
ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er; |
| 1008 |
dGk=-IC_Gk*dt/IC_tauGk +IC_b*s; |
|
| 1032 |
dGk=-IC_Gk*ANdt/IC_tauGk +IC_b*s;
|
|
| 1009 | 1033 |
IC_Gk=IC_Gk+dGk; |
| 1010 | 1034 |
|
| 1011 | 1035 |
% After a spike, the threshold is raised |
| 1012 | 1036 |
% otherwise it settles to its baseline |
| 1013 |
dTh=-(IC_Th-Th0)*dt/IC_tauTh +s*IC_c; |
|
| 1037 |
dTh=-(IC_Th-Th0)*ANdt/IC_tauTh +s*IC_c;
|
|
| 1014 | 1038 |
IC_Th=IC_Th+dTh; |
| 1015 | 1039 |
end |
| 1016 | 1040 |
end |
| ... | ... | |
| 1029 | 1053 |
lastCell=nCellsPerTau; |
| 1030 | 1054 |
for tauCount=1:nANfiberTypes |
| 1031 | 1055 |
% separate rates according to fiber types |
| 1056 |
% currently only the last segment is saved |
|
| 1032 | 1057 |
ICfiberTypeRates(tauCount, ... |
| 1033 | 1058 |
reducedSegmentPTR:shorterSegmentEndPTR)=... |
| 1034 | 1059 |
sum(ICspikes(firstCell:lastCell, :))... |
| ... | ... | |
| 1036 | 1061 |
firstCell=firstCell+nCellsPerTau; |
| 1037 | 1062 |
lastCell=lastCell+nCellsPerTau; |
| 1038 | 1063 |
end |
| 1039 |
ICoutput(:, reducedSegmentPTR:shorterSegmentEndPTR)=ICspikes; |
|
| 1040 |
|
|
| 1064 |
|
|
| 1065 |
ICoutput(:,reducedSegmentPTR:shorterSegmentEndPTR)=ICspikes; |
|
| 1066 |
|
|
| 1067 |
% store membrane output on original dt scale |
|
| 1041 | 1068 |
if nBFs==1 % single channel |
| 1042 | 1069 |
x= repmat(ICmembranePotential(1,:), ANspeedUpFactor,1); |
| 1043 | 1070 |
x= reshape(x,1,segmentLength); |
| 1044 | 1071 |
if nANfiberTypes>1 % save HSR and LSR |
| 1045 |
y= repmat(ICmembranePotential(end,:), ANspeedUpFactor,1); |
|
| 1072 |
y=repmat(ICmembranePotential(end,:),... |
|
| 1073 |
ANspeedUpFactor,1); |
|
| 1046 | 1074 |
y= reshape(y,1,segmentLength); |
| 1047 | 1075 |
x=[x; y]; |
| 1048 | 1076 |
end |
Also available in: Unified diff