% plotTraceVariableCurrentsFormal.m
% Jessica Parker, October 16, 2024
%
% This Matlab script produces Fig. 5A1, plotting measured INaP, normalized gNaP as measured, Vclamp, IPump, the actual 
% model INaPn, and the measured INaP subtracted by the model IPump. The time scale plotted goes from 10 s before the 
% step to -40 mV to 20 s after the step.

clear all;
close all;

vrsn = 'A';
leakDir = '4_0_0_1';

run1 = 4;
run2 = 0;
run3 = 0;
run4 = 2;
run5i = 1;
run5f = 3;
run6i = 1;
run6f = 3;
parname = 'I_{PumpMax}';
parunits = 'pA';
parcan = 40.0; % Ipumpmax value for run5 = 1
parstep = 40.0; % Increase in Ipumpmax between consecutive run5 instances

dataDir = 'data/';
figDir = 'plots/';
prfx = 'NaiMnap'; % Prefix used in the name of the simulation data file
nvar = 2; % Number of variables in simulation data file
Naindx = 1; % Index of [Na+]i variable in simulation data file
mNaPindx = 2; % Index of mNaP variable in simulation data file
Vref = -90; % Holding potential
tint = 0.0001; % Time step used by integrator
Vsteps = [-90, -80, -40]; % Membrane potential protocol

% Figure properties and options
fntnm = 'Lucida Sans';
fntsz = 16;
fntwght = 'bold';
vsblty = 'on'; % Set to 'off' to make figures not pop up or set to 'on'
CloseFigs = 0; % Set to 1 to automatically close figures after printing them to a file, or set to 0 to leave them open.
xmin = 70;
xmax = 100;
xl = xmax - xmin;
OrangishRed = [255, 60, 0]/255;
BrightPurple = [153 0 153]/255;
DarkTurquoise = [22, 150, 149]/255;
clrs = [DarkTurquoise; BrightPurple; OrangishRed];
SetYLimits = 1;
setymins = [-210, -0.1];
setymaxs = [0, 1.1];

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

readpars; % Model parameters read from readpars.m, not all of these are actually used by this code
assmdENa = 45; % This is assumed Na reversal potential used to measure gNaP

setyls = setymaxs-setymins;
for run5 = run5i:run5f
    bb = run5;
    Ipmpmx = parcan + (run5-1)*parstep;

    yy = [];
    leakyy = [];
    for run6=run6i:run6f   % Concatenating portions of time saved in separate files

        dir3 = [dir2 '_' num2str(run5)];

        VarN = [dataDir prfx dir3 '_' num2str(run6) '.dat'];
        fl = fopen(VarN);
        yy0 = fread(fl,[nvar, 10000000],'double')'; % Loading simulation data
        fclose(fl);

        yy = [yy(1:end-1,:); yy0];   % Remove last point here, because otherwise you will repeat the same time point
        indxs(run6) = length(yy);
        lngths(run6) = length(yy0(:,1));
        if run6 < run6f
            lngths(run6) = lngths(run6) - 1;
        end

        leakFlN = ['../' leakDir '/' dataDir prfx leakDir '_' num2str(run5) '_' num2str(run6) '.dat'];
        leakFl = fopen(leakFlN);
        leakyy0 = fread(leakFl,[nvar, 10000000],'double')'; % Loading simulation data from 4_0_0_1 used for P/4 leak subtraction
        fclose(leakFl);

        leakyy = [leakyy(1:end-1,:); leakyy0];
        leakIndxs(run6) = length(leakyy);
    end

    Vm = [];
    leakVm = [];
    for aa = run6i:run6f
        Vm = [Vm; Vsteps(aa)*ones(lngths(aa),1)]; % Clamped membrane potential over time
        leakVm = [leakVm; (Vref + 0.25*(Vsteps(aa) - Vref))*ones(lngths(aa),1)]; % Clamped membrane potential for leak measurement in 4_0_0_1
    end

    Ipump = Ipmpmx./(1.0+(Naih*exp(alpha*Vm/rtf)./yy(:,Naindx)).^nhillNa); % Ipump in 4_0_0_2
    Inap = gNaP*yy(:,mNaPindx).*(Vm-ENa); % INaP in 4_0_0_2
    IleakNa = gL*(EL-EK)/(ENa-EK)*(Vm-ENa); % ILeakNa in 4_0_0_2
    IleakK = gL*(EL-ENa)/(EK-ENa)*(Vm-EK); % ILeakK in 4_0_0_2
    Ileak = IleakNa + IleakK;
    Itot = Ipump + Inap + Ileak;
    Ipumps{bb} = Ipump; % Using cell arrays to make plotting traces together easier
    Inaps{bb} = Inap;

    IleakPump = Ipmpmx./(1.0+(Naih*exp(alpha*leakVm/rtf)./leakyy(:,Naindx)).^nhillNa); % Ipump in 4_0_0_1
    IleakNap = gNaP*leakyy(:,mNaPindx).*(leakVm-ENa); % INaP in 4_0_0_1
    IleakLNa = gL*(EL-EK)/(ENa-EK)*(leakVm-ENa); % ILeakNa in 4_0_0_1
    IleakLK = gL*(EL-ENa)/(EK-ENa)*(leakVm-EK); % ILeakK in 4_0_0_1
    IleakL = IleakLNa + IleakLK;
    IleakMsrd0 = IleakPump + IleakNap + IleakL;

    IleakRef = mean(IleakMsrd0(round(indxs(1)/2):indxs(1)-1));
    IleakMsrd = 4*(IleakMsrd0 - IleakRef) + IleakRef; % Measured leak current for P/4 leak subtraction
    IleakMsrd(1:indxs(1)) = IleakMsrd0(1:indxs(1));

    InapMsrd{bb} = Itot - IleakMsrd; % Subtracting measured leak current
    gnapMsrd{bb} = InapMsrd{bb}./(Vm - assmdENa); % Measuring gNaP using assumed Na reversal potential
    normGnap{bb} = gnapMsrd{bb}/max(gnapMsrd{bb}); % Normalized gNaP
    Nais{bb} = yy(:,Naindx); %[Na+]i, Using cell array to make plotting traces together easier
end

tt = 0.0:tint:tint*(indxs(run6f)-1);
xmini = round(xmin/tint) + 1;
xmaxi = round(xmax/tint) + 1;

f = figure('visible',vsblty);
f.PaperPositionMode = 'manual';
f.PaperUnits = 'inches';
f.Units = 'inches';
f.OuterPosition = [1.0 1.0 6.5 8.9];
f.InnerPosition = [0.25 0.25 6.0 8.4];
f.PaperPosition = [0.25 0.25 5.5 7.9];
f.RendererMode = 'manual';

titlPos = 0.12;

axes('position',[0.07 0.855 0.92 0.14]);
hold on;
text(0.2,0.4,'Model A','fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
for bb = run5i:run5f
    plot([0.515 0.545],[0.2+0.3*(run5f-bb) 0.2+0.3*(run5f-bb)],'color',clrs(bb,:),'linewidth',2.0);
    text(0.565,0.2+0.3*(run5f-bb),[parname ' = ' num2str(parcan+parstep*(bb-1)) ' ' parunits],'fontname',fntnm,'fontsize',fntsz-1,'fontweight',fntwght);
end
rectangle('position',[0.49, 0, 0.51, 1],'edgecolor','k','facecolor','none','linewidth',3);
xlim([0 1]);
ylim([0 1]);
axis off;

for bb = run5i:run5f
    ymins(bb) = min(InapMsrd{bb}(xmini:xmaxi));
    ymaxs(bb) = max(InapMsrd{bb}(xmini:xmaxi));
end
ymin = min(ymins);
ymax = max(ymaxs);
yl = ymax - ymin;
minmaxs(1,:) = [ymin, ymax];
InapMinLbl = -180;
InapMaxLbl = 0;

axes('position',[0.05 0.73 0.9 0.1]);
hold on;
for bb = run5f:-1:run5i
    plot(tt(xmini:xmaxi),InapMsrd{bb}(xmini:xmaxi),'color',clrs(bb,:),'linewidth',2.0);
end
xlim([xmin-0.1*xl xmax]);
plot([tt(indxs(2))+0.04*xl tt(indxs(2))+0.075*xl],InapMaxLbl*[1 1],'k','linewidth',3.5);
text(tt(indxs(2))+0.095*xl,InapMaxLbl,[num2str(InapMaxLbl) ' pA'],'fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
plot([tt(indxs(2))-0.075*xl tt(indxs(2))-0.04*xl],InapMinLbl*[1 1],'k','linewidth',3.5);
text(tt(indxs(2))-0.305*xl,InapMinLbl,[num2str(InapMinLbl) ' pA'],'fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
if SetYLimits
    ylim([setymins(1) setymaxs(1)]);
    text(xmin-0.9*titlPos*xl,setymins(1)+0.7*setyls(1),'I_{NaP}','fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
else
    ylim([ymin-0.1*yl ymax+0.1*yl]);
    text(xmin-0.9*titlPos*xl,ymin+0.5*yl,'I_{NaP}','fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
end
axis off;

for bb = run5i:run5f
    ymins(bb) = min(normGnap{bb}(xmini:xmaxi));
    ymaxs(bb) = max(normGnap{bb}(xmini:xmaxi));
end
ymin = min(ymins);
ymax = max(ymaxs);
yl = ymax - ymin;
minmaxs(2,:) = [ymin, ymax];

axes('position',[0.05 0.59 0.9 0.1]);
hold on;
for bb = run5f:-1:run5i
    plot(tt(xmini:xmaxi),normGnap{bb}(xmini:xmaxi),'color',clrs(bb,:),'linewidth',2.0);
end
xlim([xmin-0.1*xl xmax]);
plot([tt(indxs(2))-0.075*xl tt(indxs(2))-0.04*xl],[1 1],'k','linewidth',3.5);
text(tt(indxs(2))-0.13*xl,1,'1','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
plot([tt(indxs(2))+0.04*xl tt(indxs(2))+0.075*xl],[0 0],'k','linewidth',3.5);
text(tt(indxs(2))+0.095*xl,0,'0','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
if SetYLimits
    ylim([setymins(2) setymaxs(2)]);
    text(xmin-titlPos*xl,setymins(2)+0.6*setyls(2),'<g_{NaP}>','fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
else
    ylim([ymin-0.1*yl ymax+0.1*yl]);
    text(xmin-titlPos*xl,ymin+0.6*yl,'<g_{NaP}>','fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
end
axis off;

ymin = min(Vsteps(2:end));
ymax = max(Vsteps(2:end));
yl = ymax - ymin;
minmaxs(3,:) = [ymin, ymax];

axes('position',[0.05 0.475 0.9 0.075]);
hold on;
plot(tt(xmini:xmaxi),Vm(xmini:xmaxi),'k','linewidth',2.0);
xlim([xmin-0.1*xl xmax]);
text(xmin-titlPos*xl,ymin+0.5*yl,'V_{clamp}','fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
plot([tt(indxs(2))-0.075*xl tt(indxs(2))-0.04*xl],[-40 -40],'k','linewidth',3.5);
text(tt(indxs(2))-0.29*xl,-40,'-40 mV','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
plot([tt(indxs(2))+0.04*xl tt(indxs(2))+0.075*xl],[-80 -80],'k','linewidth',3.5);
text(tt(indxs(2))+0.095*xl,-80,'-80 mV','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
plot([xmax-8 xmax-3],[-67 -67],'k','linewidth',3.5);
text(xmax-6.6,-80,'5 s','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
ylim([ymin-0.1*yl ymax+0.1*yl]);
axis off

for bb = run5i:run5f
    ymins(bb) = min(Ipumps{bb}(xmini:xmaxi));
    ymaxs(bb) = max(Ipumps{bb}(xmini:xmaxi));
end
ymin = min(ymins);
ymax = max(ymaxs);
yl = ymax - ymin;
minmaxs(4,:) = [ymin, ymax];
IpmpMaxLbl = 60;
IpmpMinLbl = 10;

axes('position',[0.05 0.33 0.9 0.1]);
hold on;
for bb = run5f:-1:run5i
    plot(tt(xmini:xmaxi),Ipumps{bb}(xmini:xmaxi),'color',clrs(bb,:),'linewidth',2.0);
end
xlim([xmin-0.1*xl xmax]);
plot([tt(indxs(2))-0.075*xl tt(indxs(2))-0.04*xl],IpmpMaxLbl*[1 1],'k','linewidth',3.5);
text(tt(indxs(2))-0.275*xl,IpmpMaxLbl,[num2str(IpmpMaxLbl) ' pA'],'fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
plot([tt(indxs(2))+0.04*xl tt(indxs(2))+0.075*xl],IpmpMinLbl*[1 1],'k','linewidth',3.5);
text(tt(indxs(2))+0.095*xl,IpmpMinLbl,[num2str(IpmpMinLbl) ' pA'],'fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
ylim([ymin-0.1*yl ymax+0.1*yl]);
text(xmin-titlPos*xl,ymin+0.5*yl,'I_{Pump}','fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
axis off;

for bb = run5i:run5f
    ymins(bb) = min(Inaps{bb}(xmini:xmaxi));
    ymaxs(bb) = max(Inaps{bb}(xmini:xmaxi));
end
ymin = min(ymins);
ymax = max(ymaxs);
yl = ymax - ymin;
minmaxs(6,:) = [ymin, ymax];

axes('position',[0.05 0.19 0.9 0.1]);
hold on;
for bb = run5f:-1:run5i
    plot(tt(xmini:xmaxi),Inaps{bb}(xmini:xmaxi),'color',clrs(bb,:),'linewidth',2.0);
end
xlim([xmin-0.1*xl xmax]);
plot([tt(indxs(2))-0.075*xl tt(indxs(2))-0.04*xl],InapMinLbl*[1 1],'k','linewidth',3.5);
text(tt(indxs(2))-0.305*xl,InapMinLbl,[num2str(InapMinLbl) ' pA'],'fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
plot([tt(indxs(2))+0.04*xl tt(indxs(2))+0.075*xl],InapMaxLbl*[1 1],'k','linewidth',3.5);
text(tt(indxs(2))+0.095*xl,InapMaxLbl,[num2str(InapMaxLbl) ' pA'],'fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
ylim([ymin-0.1*yl ymax+0.1*yl]);
text(xmin-titlPos*xl,ymin+0.7*yl,'I_{NaPn}','fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
axis off;

Ichck = InapMsrd{bb} - Ipumps{bb};
for bb = run5i:run5f
    ymins(bb) = min(Ichck(xmini:xmaxi));
    ymaxs(bb) = max(Ichck(xmini:xmaxi));
end
ymin = min(ymins);
ymax = max(ymaxs);
yl = ymax - ymin;
minmaxs(5,:) = [ymin, ymax];
IchckMaxLbl = 0;
IchckMinLbl = -180;

axes('position',[0.05 0.05 0.9 0.1]);
hold on;
for bb = run5f:-1:run5i
    plot(tt(xmini:xmaxi),Ichck(xmini:xmaxi),'color',clrs(bb,:),'linewidth',2.0);
end
xlim([xmin-0.1*xl xmax]);
plot([tt(indxs(2))-0.075*xl tt(indxs(2))-0.04*xl],IchckMinLbl*[1 1],'k','linewidth',3.5);
text(tt(indxs(2))-0.305*xl,IchckMinLbl,[num2str(IchckMinLbl) ' pA'],'fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
plot([tt(indxs(2))+0.04*xl tt(indxs(2))+0.075*xl],IchckMaxLbl*[1 1],'k','linewidth',3.5);
text(tt(indxs(2))+0.095*xl,IchckMaxLbl,[num2str(IchckMaxLbl) ' pA'],'fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
ylim([ymin-0.1*yl ymax+0.1*yl]);
text(xmin-1.2*titlPos*xl,ymin+0.7*yl,'I_{NaP}- I_{Pump}','fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
axis off;

print(f,[figDir 'multipleVariablesCurrents' dir2 vrsn '.eps'],'-depsc','-r0');

if strcmp(vsblty,'off')
    close(f);
elseif CloseFigs
    close(f);
end