% plotReducedModelBifurcation.m
% Jessica Parker, October 16, 2024
%
% This Matlab script plots the min and max membrane potentials versus [Na+]i of the reduced constant-[Na+]i model, where
% [Na+]i has been held constant and varied as a parameter. The data of the full model is also plotted, overlayed on top
% of the bifurcation diagram of the reduced model. The phase of the full model data is color coded as shown in Vm vs time
% plotted beneath the bifurcation diagram.
clear all;
close all;
vrsn = 'J';
fntnm = 'Lucida Sans';
fntsz = 16;
fntwght = 'bold';
pgHght = 5.25;
dir = '2_0_0';
run4s = [1, 2];
PlotTrace = 1;
MarkActivityType = 0;
parcans = [13, 35]; % The starting value of [Na+]i for 2_0_0_1 and 2_0_0_2 respectively
parsteps = [0.25, -0.25]; % The step size for [Na+]i between consecutive run5 instances for 2_0_0_1 and 2_0_0_2 respectively
traceDir = ['../FullCanonical/data/']; % Directory where full model simulation data is held
traceName = '1_0_0_0_0'; % The run numbers for the full model simulation data
MarkActivityType = 1; % Set to 1 to use different marker for different activity types, like silent versus spiking
ShowLegend = 0; % Set to 1 one to include a legend for the different marker types
nvar = 5; % Number of variables in the full model simulation file (1_0_0_0_0 in /FullCanonical/)
parname = '[Na^+]_i (mM)';
ipar = 4; % The index of the variable in the full model corresponding to the parameter varied here, [Na+]i
nbcs = 3; %% Number of bursts shown on the bottom colored trace
SetXLim = 1; % Set to 1 to control the x-axis range, or set to 0 to set it automatically
xmin = 13; %% Only matters if SetXLim is set to 1
xmax = 35; %% Only matters if SetXLim is set to 1
SetYLim = 1; % Set to 1 to control the x-axis range, or set to 0 to set it automatically
yminSet = -90; %% Only matters if SetYLim is set to 1
ymaxSet = 12; %% Only matters if SetYLim is set to 1
if PlotTrace
tint = 0.0001;
VarN=[traceDir 'Vall' traceName '_1.dat']; % Name of file with full model simulation data
f1=fopen(VarN);
yy=fread(f1,[nvar,10000000],'double')'; % Load full model simulation data
tt = 0.0:tint:(length(yy(:,1))-1)*tint;
end
% Creating cell matrices that will be filled later
maxVs = cell(length(run4s));
minVs = maxVs;
ActvtyTyps = maxVs;
pars = maxVs;
% Creating arrays that will be filled later
lns = zeros(1,length(run4s));
minPars = lns;
maxPars = lns;
totminV = lns;
totmaxV = lns;
for aa=1:length(run4s)
maxVs{aa} = load([dir '_' num2str(run4s(aa)) '/data/maxVs' dir '_' num2str(run4s(aa)) 'sc.txt']); % Maximum membrane potential
minVs{aa} = load([dir '_' num2str(run4s(aa)) '/data/minVs' dir '_' num2str(run4s(aa)) 'sc.txt']); % Minimum membrane potential
ActvtyTyps{aa} = load([dir '_' num2str(run4s(aa)) '/data/type' dir '_' num2str(run4s(aa)) 'sc.txt']); % Activity type, silent, spiking,
% depolarization block, or bursting
lns(aa) = length(maxVs{aa}); % Number of run5 instances for this run4
pars{aa} = parcans(aa):parsteps(aa):parcans(aa)+(lns(aa)-1)*parsteps(aa); % Array of parameter values ([Na+]i)
minPars(aa) = min(pars{aa}); % Minimum parameter value
maxPars(aa) = max(pars{aa}); % Maximum parameter value
totminV(aa) = min(minVs{aa}); % Overall minimum of all the minimum membrane potentials
totmaxV(aa) = max(maxVs{aa}); % Overall maximum of all the maximum membrane potentials
end
minPar = min(minPars); % Overall minimum parameter value
maxPar = max(maxPars); % Overall maximum parameter value
difPar = maxPar-minPar; % Width of the range of parameter values used
% Finding min and max for the y-axis needed if the y-axis range is set automatically
ymin0 = min(totminV);
ymax0 = max(totmaxV);
if PlotTrace
ymin = min(ymin0,min(yy(:,1)));
ymax = max(ymax0,max(yy(:,1)));
else
ymin = ymin0;
ymax = ymax0;
end
ydif = ymax-ymin;
if PlotTrace
nclrs = 50; % Total number of colors used to plot the full model trace
nclrsIBI = round(0.5*nclrs); % Number of colors used to split up the IBI of the full model trace
nclrsBD = nclrs-nclrsIBI; % Number of colors used to split up the BD of the full model trace
bc = load([traceDir 'bd1j' traceName '_1sc.txt']); % burst characteristics of the full model bursting
mnBP = mean(bc(2,:)); % mean burst period of the full model
for bb = 1:nbcs+1
ttclr = tt(round(bc(1,bb)/tint):round((bc(1,bb)+bc(3,bb))/tint)); % Time points during burst duration
tt0 = tt(round((bc(1,bb)+bc(3,bb))/tint)+1:round(bc(1,bb+1)/tint)); % Time points during interburst interval
% Dulling and darkening the colors slightly so they aren't so bright.
clrsX = hsv(nclrs) - [0.2, 0.2, 0.2];
for cc = 1:3
negtv = find(clrsX(:,cc) < 0);
clrsX(negtv,cc) = 0;
end
clrs = clrsX([nclrs:-1:1],:);
% Sectioning the BD into [nclrsBD] sections.
for aa = 1:nclrsBD
iclrs{aa} = round(bc(1,bb)/tint)+(aa-1)*round(length(ttclr)/nclrsBD+1):round(bc(1,bb)/tint)+aa*round(length(ttclr)/nclrsBD+1);
end
% Sectioning the IBI into [nclrsIBI] sections.
for aa = 1:nclrsIBI
iclrs{aa+nclrsBD} = round(bc(1,bb)/tint)+nclrsBD*round(length(ttclr)/nclrsBD+1)+[(aa-1)*round(length(tt0)/nclrsIBI+1):aa*round(length(tt0)/nclrsIBI+1)];
end
ttclrs{bb} = ttclr;
clrsXs{bb} = clrsX;
allclrs{bb} = clrs;
alliclrs{bb} = iclrs;
tt0s{bb} = tt0;
end
end
f = figure();
f.PaperPositionMode = 'manual';
f.PaperUnits = 'inches';
f.Units = 'inches';
f.OuterPosition = [1 1 1.1*pgHght+1 pgHght+ShowLegend*0.25*pgHght+1];
f.InnerPosition = [0.25 0.25 1.1*pgHght+0.5 pgHght+ShowLegend*0.25*pgHght+0.5];
f.PaperPosition = [0.25 0.25 1.1*pgHght pgHght+ShowLegend*0.25*pgHght];
f.RendererMode = 'manual';
axes('position',[0.2 0.3 0.75 0.67]);
hold on;
for bb=1:length(run4s)
for aa = 1:length(pars{bb})
if ~MarkActivityType
plot(pars{bb}(aa),maxVs{bb}(aa),'k.','markersize',12);
plot(pars{bb}(aa),minVs{bb}(aa),'k.','markersize',12);
else
if ActvtyTyps{bb}(aa) == 1 % Bursting, marked by open circle
plot(pars{bb}(aa),maxVs{bb}(aa),'ko','markersize',7,'linewidth',2);
plot(pars{bb}(aa),minVs{bb}(aa),'ko','markersize',7,'linewidth',2);
elseif ActvtyTyps{bb}(aa) == 2 % Spiking, marked by dot
plot(pars{bb}(aa),maxVs{bb}(aa),'k.','markersize',12,'linewidth',2);
plot(pars{bb}(aa),minVs{bb}(aa),'k.','markersize',12,'linewidth',2);
elseif ActvtyTyps{bb}(aa) == 4 % Depolarization Block, marked by x
plot(pars{bb}(aa),maxVs{bb}(aa),'kx','markersize',7,'linewidth',2);
plot(pars{bb}(aa),minVs{bb}(aa),'kx','markersize',7,'linewidth',2);
elseif ActvtyTyps{bb}(aa) == 3 % Silence, marked by triangle
plot(pars{bb}(aa),maxVs{bb}(aa),'k^','markersize',4.5,'markerfacecolor',[0 0 0]);
plot(pars{bb}(aa),minVs{bb}(aa),'k^','markersize',4.5,'markerfacecolor',[0 0 0]);
end
p1 = plot(0,ymin-100,'k.','markersize',12);
p2 = plot(1,ymin-100,'kx','markersize',8,'linewidth',2);
p3 = plot(2,ymin-100,'ko','markersize',7,'linewidth',2);
p4 = plot(3,ymin-100,'k^','markersize',7,'markerfacecolor',[0 0 0]);
end
end
end
% Plotting full model data on top of bifurcation diagram of reduced model
if PlotTrace
iclrs = alliclrs{2};
clrs = allclrs{2};
for aa = 1:nclrs
plot(yy(iclrs{aa},ipar),yy(iclrs{aa},1),'color',[clrs(aa,:), 0.8],'linewidth',2);
end
end
if MarkActivityType && ShowLegend
legend([p1, p2, p3, p4],{'Silence','Spiking','Bursting','Depolarization Block'},'location','northoutside','fontsize',fntsz-2);
end
ylabel('Membrane Potential (mV)');
xlabel(parname);
if SetYLim
ylim([yminSet ymaxSet]);
else
ylim([ymin-0.1*ydif ymax+0.05*ydif]);
end
if SetXLim
xlim([xmin xmax]);
else
xlim([minPar-0.05*difPar maxPar+0.05*difPar]);
end
ax = gca;
ax.FontName = fntnm;
ax.FontSize = fntsz;
ax.FontWeight = fntwght;
ax.LineWidth = 3.0;
% Plotting full model membrane potential versus time underneath bifurcation diagram
if PlotTrace
axes('position',[0.2 0.03 0.75 0.1]);
hold on;
for bb = 1:nbcs+1
iclrs = alliclrs{bb};
clrs = allclrs{bb};
ttclr = ttclrs{bb};
for aa = 1:nclrs
plot(tt(iclrs{aa}),yy(iclrs{aa},1),'color',clrs(aa,:),'linewidth',2);
end
xlim([ttclrs{2}(1)-0.4*mnBP tt0s{nbcs}(end)+0.5*mnBP]);
axis off;
end
end
print(f,'-depsc',['plots/reducedBifurcation' dir '_' num2str(run4s(1)) 'x' num2str(run4s(2)) '.eps'],'-r0');