% plotInapGnapVclamp.m
% Jessica Parker, October 18, 2024
%
% This Matlab script makes the same measurements and plots as measureTailCurrent.m, but it does not save the 
% data. It is useful if you want to make changes to the figures without worrying about overwriting the data files.

clear all;
close all;

vrsn = 'A';

run1 = 5;
run2 = 0;
run3 = 0;
run4 = 2;
run5i = 1;
run5f = 3;
run6i = 1;
run6f = 3;
skpPltTrc = 1; % Number of run5 instances (+1) to skip between traces plotted. To skip none, set to 1.

dir1 = [num2str(run1) '_' num2str(run2) '_' num2str(run3)];
dir2 = [dir1 '_' num2str(run4)];
leakDir = [dir1 '_1']; % Data in 4_0_0_1 is used to measure leak current using P/4 leak subtraction

dataDir = 'data/';
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
tint = 0.0001; % Time step used by integrator

fntnm = 'Lucida Sans';
fntsz = 12;
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.
ZoomInOnStepUp = 0; % Set to 1 to zoom in on the depolarization step or set to zero and set xmin and xmax below to the axis range you want
xmin = 75;
xmax = 100;

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

yyStrt = [];

for run5 = run5i:skpPltTrc:run5f
    Ipmpmx = 40.0*run5;

    yy = yyStrt;
    leakyy = yyStrt;
    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 this time point will be repeating as the first time point of the next file
        
        indxs(run6) = length(yy);
        if run6 < run6f
            lngths(run6) = length(yy0(:,1))-1;
        else
            lngths(run6) = length(yy0(:,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 5_0_0_1 used for P/4 leak subtraction
        fclose(leakFl);

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

    Vsteps = [-90, -80, -40];
    VstepsLk = 0.25*(Vsteps - Vsteps(1)) + Vsteps(1);
    Vm = [];
    leakVm = [];
    for aa = 1:length(Vsteps)
        Vm = [Vm; Vsteps(aa)*ones(lngths(aa),1)]; % Clamped membrane potential over time
        leakVm = [leakVm; VstepsLk(aa)*ones(lngths(aa),1)]; % Clamped membrane potential for leak measurement in 5_0_0_1
    end

    Inap = gNaP*yy(:,mNaPindx).*(Vm-rtf*log(Nae./yy(:,Naindx))); % INaP in 5_0_0_2
    IleakNa = gL*(EL-EK)/(ENa-EK)*(Vm-rtf*log(Nae./yy(:,Naindx))); % ILeakNa in 5_0_0_2
    IleakK = gL*(EL-ENa)/(EK-ENa)*(Vm-EK); % ILeakK in 5_0_0_2
    Ileak = IleakNa + IleakK;
    Itot = Inap + Ileak;

    IleakNap = gNaP*leakyy(:,mNaPindx).*(leakVm-rtf*log(Nae./leakyy(:,Naindx))); % INaP in 5_0_0_1
    IleakLNa = gL*(EL-EK)/(ENa-EK)*(leakVm-rtf*log(Nae./leakyy(:,Naindx))); % ILeakNa in 5_0_0_1
    IleakLK = gL*(EL-ENa)/(EK-ENa)*(leakVm-EK); % ILeakK in 5_0_0_1
    IleakL = IleakLNa + IleakLK;
    IleakMsrd0 = IleakNap + IleakL;
    IleakMsrd = IleakMsrd0;
    IleakIntl = mean(IleakMsrd0(round(indxs(1)/2):indxs(1)));
    IleakMsrd(indxs(1)+1:indxs(end)) = IleakIntl + 4*(IleakMsrd0(indxs(1)+1:indxs(end)) - IleakIntl); % Measured leak current for P/4 subtraction

    tt = 0.0:tint:tint*(indxs(end)-1);

    InapMsrd = Itot - IleakMsrd; % Subtracting measured leak current
    gnapMsrd = InapMsrd./(Vm - assmdENa); % Measuring gNaP using assumed Na reversal potential

    if ZoomInOnStepUp
        xmin = tt(indxs(end-1)+1)-0.1;
        xmax = tt(indxs(end-1)+1) + 0.2;
        mrk = 'zm';
    else
        mrk = '';
    end

    f = figure('visible',vsblty);
    f.PaperPositionMode = 'manual';
    f.PaperUnits = 'inches';
    f.Units = 'inches';
    f.OuterPosition = [1 1 8.0 5.0];
    f.InnerPosition = [0.25 0.25 7.5 4.5];
    f.PaperPosition = [0.25 0.25 7.0 4.0];
    f.RendererMode = 'manual';

    axes('position',[0.15 0.75 0.8 0.2]);
    hold on;
    plot(tt, InapMsrd, 'k', 'linewidth', 1.5);
    ylabel('I_{NaP} (pA)');
    xticks([]);
    xlim([xmin xmax]);
    box off;
    ax = gca;
    ax.FontName = fntnm;
    ax.FontSize = fntsz;
    ax.LineWidth = 2.0;

    axes('position',[0.15 0.46 0.8 0.2]);
    hold on;
    plot(tt, gnapMsrd, 'k', 'linewidth', 1.5);
    ylabel('g_{NaP} (nS)');
    xticks([]);
    xlim([xmin xmax]);
    box off;
    ax = gca;
    ax.FontName = fntnm;
    ax.FontSize = fntsz;
    ax.LineWidth = 2.0;

    axes('position',[0.15 0.17 0.8 0.2]);
    hold on;
    plot(tt, Vm, 'k', 'linewidth', 1.5);
    ylabel('V_{clamp} (mV)');
    xlabel('Time (s)');
    xlim([xmin xmax]);
    box off;
    ax = gca;
    ax.FontName = fntnm;
    ax.FontSize = fntsz;
    ax.LineWidth = 2.0;

    print(f,['plots/InapGnapVclamp' dir2 '_' num2str(run5) vrsn mrk '.eps'],'-depsc','-r0');

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