%plotVNaiIpumpInap.m
%Jessica Parker, October 10, 2024
%
%This matlab script is used to produce a simple plot of the membrane potential, [Na+]i, Ipump, and Inap over time. 

clear all;
close all;

mrk = 'sc';

run1 = 1;
run2 = 0;
run3 = 0;
run4i = 0;
run4f = 0;
run5i = 0;
run5f = 0;
run6i = 1;
run6f = 1;
datadir = 'data/';
prfx = 'Vall';
nvar = 5;
vsblty = 'on';
CloseFigs = 1;
SingleBurst = 0;
MarkSpikeThresh = 0;
DisplayActivityType = 0;
DisplayBurstChars = 1;
figMrk = '';

xmin0 = 0; %Set the time range that you want to plot
xmax0 = 20;

%% ShowParameterAsTitle and OtherTitle cannot both be 1, but they can both be 0
ShowParameterAsTitle = 0; %boolean, set to 1 if you want a parameter value in the title
parname = 'g_{NaP}'; %Name of parameter to show in the title
parunits = 'nS'; %Units of parameter to show in the title
parcan = 3.65; %set to initial parameter value in sweep
parstep = 0.05; %set to parameter step size in sweep
titlExtra = ''; %% Anything extra that you want to put at the end of the title. Set to '' if you don't want to add anything.

Sweeping2Pars = 0; %% Use this by setting it to 1 if ShowParameterAsTitle is set to 1 and there are 2 parameters being varied.
SweepingIndependently = 0; %% Set SweepingIndependently to 1 if the two parameters are swept independently and referenced by different run numbers.
parname2 = 'I_{PumpMax}'; %% If SweepingIndependently is set to 1, the first parameter that you set as "parname" should
parunits2 = 'pA'; %% be the one represented by run4 and the second that you set as "parname2" should be the one represented by run5.
parcan2 = 50;
parstep2 = 0.5;

readpars; %% Calls to readpars.m to read the parameter values.

OtherTitle = 0; %% Use this if you want to have something else as the title where all plots will have the same title.
titlText = ''; %% Text you want to put in the title of all plots

nFigRows = 4;
titles = {'V_m (mV)','[Na^+]_i (mM)','I_{Pump} (pA)','I_{NaPn} (pA)'};

pbp = 0.1; %If you set SingleBurst to 1, this is the proportion of burst period to plot on either side of the burst 
           %duration.
tint = 0.0001; %Time step between data points
fntsz = 14; %Font size
fntnm = 'Lucida Sans';

if contains(prfx,'Vall')
  nNai = 4;
  mNaPi = 5;
  hNaFi = 3;
  mKi = 2;
elseif contains(prfx,'VNa')
  nNai = 2;
else
    disp('Check variables saved in your simulation data file.');
end

dir2 = [num2str(run1) '_' num2str(run2) '_' num2str(run3)];

aa = 1;
for run4=run4i:run4f
  typesFound = 0;
  if DisplayActivityType && exist([datadir 'type' dir2 '_' num2str(run4) mrk '.txt'])
    typ = load([datadir 'type' dir2 '_' num2str(run4) mrk '.txt']);
    typesFound = 1;
    types = {'Bursting','Spiking','Silence','Depolarization Block','Depolarized Oscillations','Unknown','Plateau Bursting'};
    mrk3 = 'Bd';
  else
    mrk3 = '';
  end

  for run5=run5i:run5f
    
    yy = [];
    for run6=run6i:run6f   %Concatenating portions of time saved in separate files
      dir3 = [dir2 '_' num2str(run4) '_' num2str(run5)];
      
      VarN=[datadir prfx dir3 '_' num2str(run6) '.dat'];
      f1=fopen(VarN);
      yy1=fread(f1,[nvar, 10000000],'double')';
      fclose(f1);
      
      yy = [yy; yy1(1:end-1,:)];   %Remove last point here, because otherwise
    end                              %you will repeat the same time point
				%twice between consecutive run6 files
    
    lnt = length(yy(:,1));
    tt = 0:tint:tint*(lnt-1);

    Ipump = Ipmpmx./(1.0+(Naih*exp(-1*alpha*yy(:,1)/rtf)./yy(:,nNai)).^nhillNa);    
    Inap =  gNaP*yy(:,mNaPi).*(yy(:,1)-rtf*log(Nae./yy(:,nNai)));
    
    plotDat = zeros(length(yy(:,1)),2);
    plotDat(:,1) = yy(:,1)';
    plotDat(:,2) = yy(:,nNai)';
    plotDat(:,3) = Ipump';
    plotDat(:,4) = Inap';
    
    bc = [];
    if SingleBurst 
      if exist([datadir 'bd1j' dir3 '_' num2str(run6f) mrk '.txt'])
        bc = load([datadir 'bd1j' dir3 '_' num2str(run6f) mrk '.txt']);
        xmin = bc(1,2)-pbp*bc(2,2);
        xmax = bc(1,2)+bc(3,2)+pbp*bc(2,2);
      elseif exist([datadir 'plts' dir3 '_' num2str(run6f) mrk '.txt'])
        bc = load([datadir 'plts' dir3 '_' num2str(run6f) mrk '.txt']);
        xmin = bc(1,2)-pbp*bc(2,2);
        xmax = bc(1,2)+bc(3,2)+pbp*bc(2,2);
      else 
        xmin = xmin0;
        xmax = xmax0;
      end
    else
      xmin = xmin0;
      xmax = xmax0;
    end

    BurstingFound = 0;
    if DisplayBurstChars 
      if exist([datadir 'bd1j' dir3 '_' num2str(run6f) mrk '.txt'])
        bc = load([datadir 'bd1j' dir3 '_' num2str(run6f) mrk '.txt']);
	BurstingFound = 1;
      elseif exist([datadir 'plts' dir3 '_' num2str(run6f) mrk '.txt'])
        bc = load([datadir 'plts' dir3 '_' num2str(run6f) mrk '.txt']);
	BurstingFound = 1;
      end
    end
    
    pgBttm = 0.75;
    if DisplayActivityType || DisplayBurstChars
      pgTop = 0.75;
    else
      pgTop = 0.5;
    end
    figSep = 0.2;
    figHght = 1.25;
    pgHght = nFigRows*figHght + pgBttm + pgTop + (nFigRows-1)*figSep;
    
    f = figure('visible',vsblty);
    f.PaperPositionMode = 'manual';
    f.PaperUnits = 'inches';
    f.Units = 'inches';
    f.OuterPosition = [1.0 1.0 8.5 pgHght+1.0];
    f.InnerPosition = [0.25 0.25 8 pgHght+0.5];
    f.PaperPosition = [0.25 0.25 7.5 pgHght];
    f.RendererMode = 'manual';
   
    figHght0 = (pgHght-pgBttm-pgTop-(nFigRows-1)*figSep)/nFigRows;

    topPerc = 1 - pgTop/pgHght;
    bttmPerc = pgBttm/pgHght;
    restPerc = topPerc - bttmPerc;
    sepPerc = figSep/pgHght;

    if DisplayActivityType && typesFound
      axes('position',[0.7 topPerc 0.3 1-topPerc]);
      text(0,0.75,types{typ(run5-run5i+1)},'fontname',fntnm,'fontsize',fntsz-3);
      axis([0 1 0 1]);
      axis off;
    elseif DisplayActivityType && ~typesFound
      disp(['Could not find activity types file for run4 = ' num2str(run4)]);
    end

    if DisplayBurstChars && BurstingFound
      axes('position',[0.2 topPerc 0.7 1-topPerc]);
      bcText = ['BP = ' num2str(round(mean(bc(2,:)),2)) ' s, BD = ' num2str(round(mean(bc(3,:)),3)) ' s, IBI = ' num2str(round(mean(bc(4,:)),2)) ' s'];
      text(0,0.6,bcText,'fontname',fntnm,'fontsize',fntsz);
      axis([0 1 0 1]);
      axis off;
    end

    maxNas(run5-run5i+1) = max(yy(:,nNai));
    minNas(run5-run5i+1) = min(yy(:,nNai));
    
    for bb = 1:nFigRows	
      axes('position',[0.15 topPerc-(restPerc/nFigRows)*bb 0.8 restPerc/nFigRows-sepPerc]);
      plot(tt,plotDat(:,bb),'k-','linewidth',1.5);
      hold on;
      if MarkSpikeThresh
        plot([xmin xmax],[-40 -40],'r--','linewidth',1);
      end
      xlim([xmin xmax]);
      if max(plotDat(:,bb)) - min(plotDat(:,bb)) > 0.01*abs(mean(plotDat(:,bb)))
	ymin = min(plotDat(:,bb));
	ymax = max(plotDat(:,bb));
	yl = ymax - ymin;
	ylim([ymin-0.1*yl ymax+0.1*yl]);
	if bb == 2
	  ylim([ymin-0.2*yl ymax+0.2*yl]);
	end
      else
	ymin = mean(plotDat(:,bb))-0.05*abs(mean(plotDat(:,bb)));
	ymax = mean(plotDat(:,bb))+0.05*abs(mean(plotDat(:,bb)));
	ylim([ymin ymax]);
      end
      ylabel(titles{bb});
      if bb == 1
	if ShowParameterAsTitle
	  if Sweeping2Pars && ~SweepingIndependently
            title([parname ' = ' num2str(parcan+parstep*(run5-1)) ' ' parunits ', ' parname2 ' = ' num2str(parcan2+parstep2*(run5-1)) ' ' parunits2 titlExtra]);
	  elseif Sweeping2Pars
	    title([parname ' = ' num2str(parcan+parstep*(run4-1)) ' ' parunits ', ' parname2 ' = ' num2str(parcan2+parstep2*(run5-1)) ' ' parunits2 titlExtra]);
	  else
	    title([parname ' = ' num2str(parcan+parstep*(run5-1)) ' ' parunits titlExtra]);
	  end
	elseif OtherTitle
	  title(titlText);
	end
      end
      if bb == nFigRows
	xlabel('Time (s)');
      else
        xticks([]);
      end
      box off;
      ax = gca;
      ax.FontName = fntnm;
      ax.FontSize = fntsz;
      ax.LineWidth = 2.0;
    end
    
    if SingleBurst && length(bc)>0
      print(f,['plots/sbVNaiIpumpInap' dir3 '_' num2str(run6) mrk3 figMrk '.eps'],'-depsc','-r0');
    else
      print(f,['plots/VNaiIpumpInap' dir3 '_' num2str(run6) mrk3 figMrk '.eps'],'-depsc','-r0');
    end
    
    if strcmp(vsblty,'off')
      close(f);
    elseif CloseFigs
      close(f);
    end
    aa = aa + 1;
  end
end