%runburstanalysisSC.m
%Jessica Parker, October 10, 2024
%
%Matlab script for calculating and saving burst characteristics. Calls to burstanalysis() function for the actual
%measurements. Plots activity with characteristics marked if plt = 1.
clear all;
mrk = 'sc';
run1 = 2;
run2 = 0;
run3i = 0;
run3f = 0;
run4i = 2; %Used as beginning of for loop to run through all run4 values
run4f = 2; %Used as end of for loop, set to last run4 value
run5i = 1; %Used as beggining of for loop
run5f = 89; %Used as end of for loop, set to last run6 value
run6i = 1;
run6f = 1;
datadir = 'data/'; %Directory where the simulation data file is held.
prfx = 'V'; %Prefix of the simulation data file name.
nvar = 1;
naIndx = 0; %Set to index of [Na+]i data in simulation file or set to zero if you don't have it or dont want to measure mean [Na+]i
AddDataToEndOfFile = 0;
plt = 0; %Set to 1 if you want to plot important aspects of episode analysis in order to check that they were calculated
%correctly. Set to 0 if you don't want to make any figures.
PlotBurstingOnly = 1;
vsblty = 'on';
CloseFigs = 1;
figDir = 'plots/'; %If you set plt to 1, this is where the figures will be saved.
ibith = 0.3; %Threshold for interburst interval, use 0.3 for full model, use 0.1 or 0.12 for Iinj-induced bursting
spth = -40.0; %Spike threshold
xmin = 0.0; %minimum time value plotted
xmax = 10.0; %maximum time value plotted
tint = 0.0001; %time step between data points
writeDirFiles = 1; %Set to 1 if you want to save files with lists of various means for dir2 or [run1]_[run2]_[run3]_[run4],
%see bottom of this file
for run3=run3i:run3f
dir = [num2str(run1) '_' num2str(run2) '_' num2str(run3)];
bcs = [];
bcsPlts = [];
for run4=run4i:run4f
dir2 = [dir '_' num2str(run4)];
cc = run4-run4i+1;
bb = 0;
dd = 0;
if naIndx > 0
mnNais = [];
minNais = [];
maxNais = [];
end
mnVs = [];
bcs = [];
bcsPlts = [];
actvtyTyp = zeros(1,run5f-run5i+1);
for run5=run5i:run5f
dir3 = [dir2 '_' num2str(run5)];
aa = run5-run5i+1;
V1 = [];
Nai = [];
for run6 = run6i:run6f %Concatenate all time portions of a simulation
dir4 = [dir3 '_' num2str(run6)];
VarN=[datadir prfx dir4 '.dat'];
f1=fopen(VarN);
yy1=fread(f1,[nvar,10000000],'double')'; %Use
%above array is [number of variables saved, any large number
%greater than the number of data points]
fclose(f1);
V1 = [V1; yy1(1:end-1,1)];
if naIndx > 0
Nai = [Nai; yy1(1:end-1,naIndx)];
end
end
lnt = length(V1);
tt = 0:tint:tint*(lnt-1);
%measure episode characteristics and burst characteristics
[bc1, tspks1, lowAmpOsc] = burstanalysis(V1,tint,spth,ibith);
if ~isempty(tspks1)
dlmwrite([datadir 'sps' dir4 mrk '.txt'], tspks1);
end
if naIndx > 0
mnNais(aa) = mean(Nai);
minNais(aa) = min(Nai);
maxNais(aa) = max(Nai);
end
mnVs(aa) = mean(V1);
maxVs(aa) = max(V1);
minVs(aa) = min(V1);
if ~isempty(tspks1)
isi = tspks1(2:end)-tspks1(1:end-1);
else
isi = [];
end
if ~isempty(bc1) && length(bc1(1,:)) > 1 %true if case is continuous bursting
if std(bc1(3,:))/mean(bc1(3,:)) > 0.2 && (std(bc1(6,:)/bc1(6,:)) > 0.2 || (mean(bc1(5,:)) < 2 && std(bc1(4,:))/mean(bc1(4,:)) > 0.1))
disp(['Chaotic spiking at run4 = ' num2str(run4) ' and run5 = ' num2str(run5)]);
actvtyTyp(aa) = 6;
elseif lowAmpOsc > 0
disp(['Plateau bursting at run4 = ' num2str(run4) ' and run5 = ' num2str(run5)]);
dlmwrite([datadir 'plts' dir4 mrk '.txt'], bc1);
actvtyTyp(aa) = 7;
bcsPlts(1,aa) = mean(bc1(2,:));
bcsPlts(2,aa) = mean(bc1(3,:));
bcsPlts(3,aa) = mean(bc1(4,:));
bcsPlts(4,aa) = mean(bc1(5,:));
bcsPlts(5,aa) = mean(bc1(6,:));
dd = dd+1;
else %Regular bursting
disp(['Bursting at run4 = ' num2str(run4) ' and run5 = ' num2str(run5)]);
% xmax = xmin+10*mean(bc1(2,:));
dlmwrite([datadir 'bd1j' dir4 mrk '.txt'], bc1);
actvtyTyp(aa) = 1;
bcs(1,aa) = mean(bc1(2,:));
bcs(2,aa) = mean(bc1(3,:));
bcs(3,aa) = mean(bc1(4,:));
bcs(4,aa) = mean(bc1(5,:));
bcs(5,aa) = mean(bc1(6,:));
bb = bb+1;
bcsPlts(1,aa) = mean(bc1(2,:));
bcsPlts(2,aa) = mean(bc1(3,:));
bcsPlts(3,aa) = mean(bc1(4,:));
bcsPlts(4,aa) = mean(bc1(5,:));
bcsPlts(5,aa) = mean(bc1(6,:));
dd = dd+1;
if plt
f = figure('visible',vsblty);
f.PaperPositionMode = 'auto';
f.PaperUnits = 'inches';
f.PaperPosition = [0 0 8 3];
if min(bc1(5,:)) > 1
xmin = bc1(1,2)-0.8*bc1(3,2);
xmax = bc1(1,2)+1.8*bc1(3,2);
else
xmin = 0;
xmax = 3*mean(bc1(2,:));
end
%xmin = 0;
%xmax = 5*mean(bc1(2,:));
%% red dots mark spikes, and cyan x's mark beginnings and ends of
%% burst durations
axes('position',[0.15 0.15 0.8 0.8]);
plot(tt,V1,'k','linewidth',1.5);
hold on;
plot(tspks1,V1(round(tspks1./tint+1)),'r.','markersize',5);
plot(bc1(1,:),-20*ones(1,length(bc1(1,:))),'bx','markersize',10,'linewidth',2);
plot(bc1(1,:)+bc1(3,:),-20*ones(1,length(bc1(1,:))),'gx','markersize',10,'linewidth',2);
ylim([-80 10]);
xlim([xmin xmax]);
ylabel('V_1 (mV)');
print(f,[figDir 'chckBC' dir3 '.png'],'-dpng','-r0');
if strcmp(vsblty,'off')
close(f);
elseif CloseFigs
close(f);
end
end
end
elseif length(tspks1) > 4 && std(isi)/mean(isi) < 0.2 %true if case is tonic spiking
%%elseif length(tspks1) > 4
disp(['Spiking but no bursting at run4 = ' num2str(run4) ' and run5 = ' num2str(run5)]);
actvtyTyp(aa) = 2;
% if length(tspks1) > 10
% xmax = xmin+tspks1(10);
% else
% xmax = run6f*100.0;
% end
if plt && ~PlotBurstingOnly
f = figure('visible',vsblty);
f.PaperPositionMode = 'auto';
f.PaperUnits = 'inches';
f.PaperPosition = [0 0 8 3];
axes('position',[0.15 0.15 0.8 0.8]);
plot(tt,V1,'k','linewidth',1.5);
hold on;
plot(tspks1,V1(round(tspks1./tint+1)),'r.','markersize',5);
ylim([-100 50]);
xlim([xmin xmax]);
ylabel('V_1 (mV)');
print(f,[figDir 'chckBC' dir3 '.png'],'-dpng','-r0');
if strcmp(vsblty,'off')
close(f);
elseif CloseFigs
close(f);
end
end
elseif length(tspks1) > 0
disp(['Chaotic spiking at run4 = ' num2str(run4) ' and run5 = ' num2str(run5)]);
actvtyTyp(aa) = 6;
if plt && ~PlotBurstingOnly
f = figure('visible',vsblty);
f.PaperPositionMode = 'auto';
f.PaperUnits = 'inches';
f.PaperPosition = [0 0 8 3];
axes('position',[0.15 0.15 0.8 0.8]);
plot(tt,V1,'k','linewidth',1.5);
hold on;
plot(tspks1,V1(round(tspks1./tint+1)),'r.','markersize',5);
ylim([-100 50]);
xlim([xmin xmax]);
ylabel('V_1 (mV)');
print(f,[figDir 'chckBC' dir3 '.png'],'-dpng','-r0');
if strcmp(vsblty,'off')
close(f);
elseif CloseFigs
close(f);
end
end
else %true if case is silent
if max(V1) < -50
disp(['Silence at run4 = ' num2str(run4) ' and run5 = ' num2str(run5)]);
actvtyTyp(aa) = 3;
else
if max(V1) - min(V1) > 5
actvtyTyp(aa) = 5;
disp(['Depolarization block oscillations at run4 = ' num2str(run4) ' and run5 = ' num2str(run5)]);
else
actvtyTyp(aa) = 4;
disp(['Depolarization block at run4 = ' num2str(run4) ' and run5 = ' num2str(run5)]);
end
end
end
end
nbrstng(cc) = bb;
if writeDirFiles
if AddDataToEndOfFile && run5i > 1
if bb>0
dlmwrite([datadir 'bcmns' dir2 mrk '.txt'],bcs,'-append');
end
if dd>0
dlmwrite([datadir 'bcPlatMns' dir2 mrk '.txt'],bcsPlts,'-append');
end
if naIndx > 0
dlmwrite([datadir 'mnNai' dir2 mrk '.txt'],mnNais','-append');
dlmwrite([datadir 'minNai' dir2 mrk '.txt'],minNais','-append');
dlmwrite([datadir 'maxNai' dir2 mrk '.txt'],maxNais','-append');
end
dlmwrite([datadir 'mnVs' dir2 mrk '.txt'],mnVs','-append');
dlmwrite([datadir 'maxVs' dir2 mrk '.txt'],maxVs','-append');
dlmwrite([datadir 'minVs' dir2 mrk '.txt'],minVs','-append');
dlmwrite([datadir 'type' dir2 mrk '.txt'],actvtyTyp','-append');
else
if bb>0
dlmwrite([datadir 'bcmns' dir2 mrk '.txt'],bcs);
end
if dd>0
dlmwrite([datadir 'bcPlatMns' dir2 mrk '.txt'],bcsPlts);
end
if naIndx > 0
dlmwrite([datadir 'mnNai' dir2 mrk '.txt'],mnNais');
dlmwrite([datadir 'minNai' dir2 mrk '.txt'],minNais');
dlmwrite([datadir 'maxNai' dir2 mrk '.txt'],maxNais');
end
dlmwrite([datadir 'mnVs' dir2 mrk '.txt'],mnVs');
dlmwrite([datadir 'maxVs' dir2 mrk '.txt'],maxVs');
dlmwrite([datadir 'minVs' dir2 mrk '.txt'],minVs');
dlmwrite([datadir 'type' dir2 mrk '.txt'],actvtyTyp');
end
end
end
end