% measureTailCurrent.m
% Jessica Parker, October 16, 2024
%
% This Matlab script measures and saves the tail current and tail conductance from the peak Gnap onwards during the last 
% membrane potential step to -40 mV.

clear all;
close all;

vrsn = 'A';
leakDir = '5_0_0_1';

run1 = 5;
run2 = 0;
run3 = 0;
run4 = 2;
run5i = 1;
run5f = 3;
run6i = 1;
run6f = 3;

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;

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

for run5 = run5i:run5f
    Ipmpmx = 40.0*run5; % Ipumpmax changes with run5

    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];
    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 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;

    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 = Itot - IleakMsrd; % Subtracting measured leak current
    gnapMsrd = InapMsrd./(Vm - assmdENa); % Measuring gNaP using assumed Na reversal potential
    mxmlGnap = max(gnapMsrd); % Maximum gNaP reached here
    normGnapMsrd = gnapMsrd/mxmlGnap; % Normalized gNaP

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

    [maxGnap, maxGindx] = max(gnapMsrd); % peak gNaP, which occurs near beginning of last step to -40 mV
    tailcurrent = InapMsrd(maxGindx:end)'; % measured INaP throughout tail (from peak gNaP onwards)
    tailcondctnc = gnapMsrd(maxGindx:end)'; % measured gNaP throughout tail
    tailnormcond = normGnapMsrd(maxGindx:end)'; % normalized gNaP throughout tail
    ttail = tt(maxGindx:end)-tt(maxGindx); % time points included in tail, reset so t=0 is the moment of peak gNaP

    dlmwrite([dataDir 'tailCurrent' dir3 vrsn '.txt'],[tailcurrent; ttail]);
    dlmwrite([dataDir 'tailConductance' dir3 vrsn '.txt'],[tailcondctnc; ttail]);
    dlmwrite([dataDir 'tailNormalizedConductance' dir3 vrsn '.txt'],[tailnormcond; ttail]);
    dlmwrite([dataDir 'tailNaConcentration' dir3 vrsn '.txt'],[yy(maxGindx:end,Naindx)'; ttail]);

    f = figure('visible',vsblty);
    f.PaperPositionMode = 'manual';
    f.PaperUnits = 'inches';
    f.Units = 'inches';
    f.OuterPosition = [1.0 1.0 7.25 4.0];
    f.InnerPosition = [0.25 0.25 6.75 3.5];
    f.PaperPosition = [0.25 0.25 6.25 3.0];
    f.RendererMode = 'manual';

    ymin = min(InapMsrd(xmini:xmaxi));
    ymax = max(InapMsrd(xmini:xmaxi));
    minmaxs(1,:) = [ymin, ymax];
    yl = ymax - ymin;

    axes('position',[0.05 0.73 0.9 0.22]);
    hold on;
    plot(tt(xmini:xmaxi),InapMsrd(xmini:xmaxi),'k','linewidth',2.0);
    xlim([xmin-0.1*xl xmax]);
    text(xmin-0.1*xl,ymin+0.6*yl,'I_{NaP}','fontname',fntnm,'fontsize',fntsz+1,'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 pA','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
    plot([tt(indxs(2))-0.075*xl tt(indxs(2))-0.04*xl],[-180 -180],'k','linewidth',3.5);
    text(tt(indxs(2))-0.28*xl,-180,'-180 pA','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
    ylim([ymin-0.1*yl ymax+0.1*yl]);
    axis off;

    ymin = min(gnapMsrd(xmini:xmaxi));
    ymax = max(gnapMsrd(xmini:xmaxi));
    minmaxs(2,:) = [ymin, ymax];
    yl = ymax - ymin;

    axes('position',[0.05 0.4 0.9 0.22]);
    hold on;
    plot(tt(xmini:xmaxi),gnapMsrd(xmini:xmaxi),'k','linewidth',2.0);
    xlim([xmin-0.1*xl xmax]);
    text(xmin-0.1*xl,ymin+0.6*yl,'g_{NaP}','fontname',fntnm,'fontsize',fntsz+1,'fontweight',fntwght);
    plot([tt(indxs(2))-0.075*xl tt(indxs(2))-0.04*xl],[1.5 1.5],'k','linewidth',3.5);
    text(tt(indxs(2))-0.245*xl,1.5,'1.5 nS','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 nS','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
    ylim([ymin-0.1*yl ymax+0.1*yl]);
    axis off;

    ymin = min(Vsteps(2:end));
    ymax = max(Vsteps(2:end));
    yl = ymax - ymin;

    axes('position',[0.05 0.07 0.9 0.22]);
    hold on;
    plot(tt(xmini:xmaxi),Vm(xmini:xmaxi),'k','linewidth',2.0);
    xlim([xmin-0.1*xl xmax]);
    text(xmin-0.1*xl,ymin+0.6*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.265*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-10 xmax-5],[-70 -70],'k','linewidth',3.5);
    text(xmax-8.5,-80,'5 s','fontname',fntnm,'fontsize',fntsz,'fontweight',fntwght);
    ylim([ymin-0.1*yl ymax+0.1*yl]);
    axis off

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

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