% clear all
% close all

% sample interval for simulation (sec)
dt = 0.1;

% whether or not to optimize the model to fit the reference data
% !!! this was NOT used for the paper
dofit = false;

% for scaling displayed reference currents
ymin = -0.08;

% model parameters (rate constants (/s or /M/s for binding rates) and interaction eneregies (kcal/mol))
ko = 0.02; % pore opening
kc = 25; % pore closing
konH = 5e4; % proton binding rate to gating site
koffH = 1; % proton unbinding rate from gating site
konB = 5e4; % Zn binding rate to blocking site
koffB = 20; % Zn unbinding rate from blocking site
konP = 35; % Zn binding rate to potentiating site
koffP = 1; % Zn unbinding rate from potentiating site
H1G1 = -4;
H1Gt = -4;
HtG1 = -1;
P1G1 = -6; % Interaction energy for open and potentiated states (kcal/mol)
P1Gt = -3; % Interaction energy for pore gate transition state in potentiated states (kcal/mol)
PtG1 = -5; % Interaction energy for potentiated site transition state in open states (kcal/mol)
H1P1 = 6;
H1Pt = 0;
HtP1 = 0;

p = [ko, kc, konH, koffH, konB, koffB, konP, koffP, H1G1, H1Gt, HtG1, P1G1, P1Gt, PtG1, H1P1, H1Pt, HtP1];

%% Fig 1 experimental data

% load experimental data
load('fig1_data_set2.mat');

% +1 so time zero starts 1 sec before pH 5.5 pulse
fig1_ref_time_s = fig1.time_ms .* 1e-3 + 1;

% shift to visually align with current responses in simulation
fig1_ref_time_s = fig1_ref_time_s - 0.1;

% assume 80 mV driving force and Po-max ~? with normalized conductance of 1 pS
fig1_Pomax = 0.45;
fig1_ref_current_pA = fig1.current_pA ./ max(max(abs(fig1.current_pA))) .* 0.08 .* fig1_Pomax;

% clip to timing of simulation
idx = (fig1_ref_time_s > -200e-6) & (fig1_ref_time_s < 12 + 200e-6);
fig1_ref_time_s = fig1_ref_time_s(idx);
fig1_ref_time_s = fig1_ref_time_s - fig1_ref_time_s(1);
fig1_ref_current_pA = fig1_ref_current_pA(idx,:);

% decimate fs=1/(200us) to simulation sample rate fs=1/dt
r = dt / 200e-6;
fig1_ref_time_ds_s = fig1_ref_time_s(1:r:end);
clear fig1_ref_current_ds_pA;
for col = 1:size(fig1_ref_current_pA, 2)
    fig1_ref_current_ds_pA(:,col) = decimate(fig1_ref_current_pA(:,col), r);
end

% clear fig1_ref_time_s fig1_ref_current_pA;

%% Fig 2 experimental data

% load experimental data
load('fig2_data.mat');

% +1 so time zero starts 1 sec before pH 5.5 pulse
fig2_refA_time_s = fig2.A.time_ms .* 1e-3 + 1;
fig2_refB_time_s = fig2.B.time_ms .* 1e-3 + 1;
fig2_refC_time_s = fig2.C.time_ms .* 1e-3 + 1;

% shift to visually align with current responses in simulation
fig2_refA_time_s = fig2_refA_time_s - 0.38;
fig2_refB_time_s = fig2_refB_time_s - 0.48;
fig2_refC_time_s = fig2_refC_time_s - 0.36;

% assume 80 mV driving force and Po-max ~? with normalized conductance of 1 pS
fig2_PomaxA = 0.1;
fig2_PomaxB = 0.5;
fig2_PomaxC = 0.9;
fig2_refA_current_pA = fig2.A.current_pA ./ max(max(abs(fig2.A.current_pA))) .* 0.08 .* fig2_PomaxA;
fig2_refB_current_pA = fig2.B.current_pA ./ max(max(abs(fig2.B.current_pA))) .* 0.08 .* fig2_PomaxB;
fig2_refC_current_pA = fig2.C.current_pA ./ max(max(abs(fig2.C.current_pA))) .* 0.08 .* fig2_PomaxC;

% clip to timing of simulation
idxA = (fig2_refA_time_s > -100e-6) & (fig2_refA_time_s < 10 + 100e-6);
idxB = (fig2_refB_time_s > -200e-6) & (fig2_refB_time_s < 10 + 200e-6);
idxC = (fig2_refC_time_s > -100e-6) & (fig2_refC_time_s < 10 + 100e-6);
fig2_refA_time_s = fig2_refA_time_s(idxA);
fig2_refB_time_s = fig2_refB_time_s(idxB);
fig2_refC_time_s = fig2_refC_time_s(idxC);
fig2_refA_time_s = fig2_refA_time_s - fig2_refA_time_s(1);
fig2_refB_time_s = fig2_refB_time_s - fig2_refB_time_s(1);
fig2_refC_time_s = fig2_refC_time_s - fig2_refC_time_s(1);
fig2_refA_current_pA = fig2_refA_current_pA(idxA,:);
fig2_refB_current_pA = fig2_refB_current_pA(idxB,:);
fig2_refC_current_pA = fig2_refC_current_pA(idxC,:);

% decimate fs=1/(200us) to simulation sample rate fs=1/dt
r = dt / 200e-6;
fig2_refA_time_ds_s = fig2_refA_time_s(1:r:end);
fig2_refB_time_ds_s = fig2_refB_time_s(1:r:end);
fig2_refC_time_ds_s = fig2_refC_time_s(1:r:end);
clear fig2_refA_current_ds_pA;
for col = 1:size(fig2_refA_current_pA, 2)
    fig2_refA_current_ds_pA(:,col) = decimate(fig2_refA_current_pA(:,col), r);
end
clear fig2_refB_current_ds_pA;
for col = 1:size(fig2_refB_current_pA, 2)
    fig2_refB_current_ds_pA(:,col) = decimate(fig2_refB_current_pA(:,col), r);
end
clear fig2_refC_current_ds_pA;
for col = 1:size(fig2_refC_current_pA, 2)
    fig2_refC_current_ds_pA(:,col) = decimate(fig2_refC_current_pA(:,col), r);
end

% clear fig2_refA_current_pA fig2_refB_current_pA fig2_refC_current_pA;

%% Fig 1 simulation

Zn = [0,0.1,0.3,1,3,10] .* 1e-3; % M
[time, I, Po, PG, PH, PB, PP, idx] = fig1_simulation(p, dt, Zn);

try
    figure(Fig1C_current);
catch
    Fig1C_current = figure;
end
hold off;
plot(fig1_ref_time_s, fig1_ref_current_pA, 'color', [.5 .5 .5], 'HandleVisibility', 'off');
hold on;
cmap = colormap(lines(7));
plot(time, I(:,1), 'color', [0 0 0], 'linewidth', 2);
set(gca, 'ColorOrder', cmap);
set(gca, 'ColorOrderIndex', 1);
plot(time, I(:,2:end), 'linewidth', 2);
xlabel('Time(s)');
ylabel('Current(pA)');
title('Fig1C: Zn dose-response');
legend({'0 mM Zn', '0.1 mM Zn', '0.3 mM Zn', '1 mM Zn', '3 mM Zn', '10 mM Zn'}, 'location', 'southwest');
if ymin
    ylim([ymin, 0.01]);
end

try
    figure(Fig1C_proba);
catch
    Fig1C_proba = figure;
end
lw = 2;
hold off;
hPo = plot(time, Po, 'b', 'linewidth', lw);
hold on;
hPG = plot(time, PG, 'c', 'linewidth', lw);
hPH = plot(time, PH, 'g', 'linewidth', lw);
% hPB = plot(time, PB, 'r', 'linewidth', lw);
hPP = plot(time, PP, 'm', 'linewidth', lw);
xlabel('Time(s)');
ylabel('Probability');
title('Fig1C: Zn dose-response');
legend([hPo(1), hPG(1), hPH(1), hPP(1)], {'Po', 'Gate', 'H site', 'P site'});

try
    figure(Fig1D);
catch
    Fig1D = figure;
end
subplot(1, 2, 1);
plot(Zn .* 1e3, min(I,[],1) ./ min(I(:,1)), 'o', 'LineWidth', 2);
set(gca, 'XScale', 'log');
xlim([1e-2, 1e2]);
xlabel('Zn (mM)');
ylabel('Potentition (fold)');
title('Fig 2D');

subplot(1, 2, 2);
preblock = I(idx{3}(1)-1,:);
postblock = max(I(idx{3},:),[],1);
plot(Zn .* 1e3, (preblock - postblock) ./ preblock, 'o', 'LineWidth', 2);
set(gca, 'XScale', 'log');
xlim([1e-2, 1e2]);
xlabel('Zn (mM)');
ylabel('Blocking (ratio)');
title('Fig 2D');

%% Fig 2 simulation

Zn = [0.3e-3, 1e-3, 3e-3]; % [M]
prepulseSec = [0, 1, 2, 4, 8, 16, 32, 64]; % Zn prepulse duration (s)
for k = 1:length(Zn)
    [time, I, Po, PG, PH, PB, PP] = fig2_simulation(p, dt, Zn(k), prepulseSec);

    idx = (time >= 65) & (time <= 75);
    time = time(idx);
    time = time - time(1);
    Po = Po(idx,:);
    PG = PG(idx,:);
    PH = PH(idx,:);
    PB = PB(idx,:);
    PP = PP(idx,:);
    I = I(idx,:);

    try
        figure(Fig2B_current(k));
    catch
        Fig2B_current(k) = figure;
    end
    hold off;
    if k == 1
        plot(fig2_refA_time_s, fig2_refA_current_pA, 'color', [.5 .5 .5], 'HandleVisibility', 'off');
    elseif k == 2
        plot(fig2_refB_time_s, fig2_refB_current_pA, 'color', [.5 .5 .5], 'HandleVisibility', 'off');
    elseif k == 3
        plot(fig2_refC_time_s, fig2_refC_current_pA, 'color', [.5 .5 .5], 'HandleVisibility', 'off');
    end
    hold on;
    plot(time, I(:,1), 'color', [0 0 0], 'linewidth', 2);
    set(gca, 'ColorOrder', cmap);
    set(gca, 'ColorOrderIndex', 1);
    plot(time, I(:,2:end), 'linewidth', 2);
    xlabel('Time(s)');
    ylabel('Current(pA)');
    title([num2str(Zn(k)*1e3) ' mM Zn Preincubation']);
    legend({'0 s', '1 s', '2 s', '4 s', '8 s', '16 s', '32 s', '64 s'}, 'location', 'southeast');
    if ymin
        ylim([ymin, 0.01]);
    end

    try
        figure(Fig2B_proba(k));
    catch
        Fig2B_proba(k) = figure;
    end
    lw = 2;
    hold off;
    hPo = plot(time, Po, 'b', 'linewidth', lw);
    hold on;
    hPG = plot(time, PG, 'c', 'linewidth', lw);
    hPH = plot(time, PH, 'g', 'linewidth', lw);
    % hPB = plot(time, PB, 'r', 'linewidth', lw);
    hPP = plot(time, PP, 'm', 'linewidth', lw);
    xlabel('Time(s)');
    ylabel('Probability');
    title([num2str(Zn(k)*1e3) ' mM Zn Preincubation']);
    legend([hPo(1), hPG(1), hPH(1), hPP(1)], {'Po', 'Gate', 'H site', 'P site'});

    try
        figure(Fig2C);
    catch
        Fig2C = figure;
    end
    if k == 1
        hold off;
    else
        hold on;
    end
    plot(1:length(prepulseSec), min(I,[],1) ./ min(I(:,1)), 'o-', 'LineWidth', 2);
    
    try
        figure(Fig2D);
    catch
        Fig2D = figure;
    end
    if k == 1
        hold off;
    else
        hold on;
    end
    [peaks, peaksIdx] = min(I,[],1);
    latencies = time(peaksIdx) - 1;
    plot(1:length(prepulseSec), latencies, 'o-', 'LineWidth', 2);
end

figure(Fig2C);
xticks(1:length(prepulseSec));
xticklabels(arrayfun(@(sec) num2str(sec), prepulseSec, 'UniformOutput', false));
xlabel('Zn exposure time (s)');
ylabel('Normalized I (I/I_{t=0})');
legend({'0.3 mM', '1 mM', '3 mM'}, 'Location', 'northoutside');

figure(Fig2D);
xticks(1:length(prepulseSec));
xticklabels(arrayfun(@(sec) num2str(sec), prepulseSec, 'UniformOutput', false));
xlabel('Zn exposure time (s)');
ylabel('Time to I_{max} (s)');
legend({'0.3 mM', '1 mM', '3 mM'}, 'Location', 'northoutside');

%% Supp Fig 1 simulation

Zn = 1e-3;
washoutSec = [0, 1, 2, 4, 8, 16, 32, 64];
[time, I, Po, PG, PH, PB, PP, idx] = supfig1_simulation(p, dt, Zn, washoutSec);

try
    figure(FigS1A_current);
catch
    FigS1A_current = figure;
end
hold off;
set(gca, 'ColorOrder', cmap);
set(gca, 'ColorOrderIndex', 1);
plot(time, I, 'linewidth', 2);
xlim([33, 43]+64-16);
xlabel('Time(s)');
ylabel('Current(pA)');
title('1 mM Zn Washout');
labels = strsplit(num2str(washoutSec));
for i = 1:length(labels)
    labels{i} = [labels{i} ' s'];
end
legend(labels, 'location', 'southwest');
if ymin
    ylim([ymin, 0.01]);
end

try
    figure(FigS1A_proba);
catch
    FigS1A_proba = figure;
end
lw = 2;
hold off;
hPo = plot(time, Po, 'b', 'linewidth', lw);
hold on;
hPG = plot(time, PG, 'c', 'linewidth', lw);
hPH = plot(time, PH, 'g', 'linewidth', lw);
% hPB = plot(time, PB, 'r', 'linewidth', lw);
hPP = plot(time, PP, 'm', 'linewidth', lw);
% xlim([33, 43]+64-16);
xlabel('Time(s)');
ylabel('Probability');
title('1 mM Zn Washout');
legend([hPo(1), hPG(1), hPH(1), hPP(1)], {'Po', 'Gate', 'H site', 'P site'});

try
    figure(FigS1BC);
catch
    FigS1BC = figure;
end
subplot(1,2,1);
peaks = min(I,[],1);
plot(1:length(washoutSec), peaks ./ min(I(:,end)), 'o-', 'LineWidth', 2);
xticks(1:length(washoutSec));
xticklabels(arrayfun(@(sec) num2str(sec), washoutSec, 'UniformOutput', false));
xlabel('Wash-off time (s)');
ylabel('Normalized I (I/I_{t=0})');
legend({'1 mM'}, 'Location', 'northoutside');

subplot(1,2,2);
[peaks, peakIdxs] = min(I,[],1);
time0 = 2+max(washoutSec)+16;
latencies = time(peakIdxs) - time0;
plot(1:length(washoutSec), latencies, 'o-', 'LineWidth', 2);
xticks(1:length(washoutSec));
xticklabels(arrayfun(@(sec) num2str(sec), washoutSec, 'UniformOutput', false));
xlabel('Wash-off time (s)');
ylabel('Time to I_{max} (s)');
legend({'1 mM'}, 'Location', 'northoutside');

%% tile windows on screen

MP = get(0, "MonitorPositions");
n = size(MP, 1); % choose last monitor
sx = MP(n,1);
sy = MP(n,2);
sw = MP(n,3);
sh = MP(n,4);
figs = [ ...
    Fig1C_current,  Fig2B_current,                  FigS1A_current    ; ...
    Fig1C_proba,    Fig2B_proba,                    FigS1A_proba      ; ...
    Fig1D,          Fig2C, Fig2D, gobjects(1),      FigS1BC    ];
nrows = size(figs,1);
ncols = size(figs,2);

% right side of screen only?
% sx = sx + sw/2;
% sw = sw/2;

w = sw / ncols;
h = sh / nrows;
for row = 1:nrows
    for col = 1:ncols
        fig = figs(row,col);
        if isgraphics(fig)
            fig.OuterPosition = [sx+(col-1)*w, sy+sh-row*h, w, h];
        end
    end
end

%% save figures

% for ext = ["svg", "png"]
% 
%     saveas(Fig1C_current, ['figs/Fig1C_current.' char(ext)]);
%     saveas(Fig1C_proba, ['figs/Fig1C_proba.' char(ext)]);
%     
%     saveas(Fig1D, ['figs/Fig1D.' char(ext)]);
%     
%     saveas(Fig2B_current(1), ['figs/Fig2B_current_300nM.' char(ext)]);
%     saveas(Fig2B_proba(1), ['figs/Fig2B_proba_300nM.' char(ext)]);
%     saveas(Fig2B_current(2), ['figs/Fig2B_current_1mM.' char(ext)]);
%     saveas(Fig2B_proba(2), ['figs/Fig2B_proba_1mM.' char(ext)]);
%     saveas(Fig2B_current(3), ['figs/Fig2B_current_3mM.' char(ext)]);
%     saveas(Fig2B_proba(3), ['figs/Fig2B_proba_3mM.' char(ext)]);
%     
%     saveas(Fig2C, ['figs/Fig2C.' char(ext)]);
%     saveas(Fig2D, ['figs/Fig2D.' char(ext)]);
%     
%     saveas(FigS1A_current, ['figs/FigS1A_current.' char(ext)]);
%     saveas(FigS1A_proba, ['figs/FigS1A_proba.' char(ext)]);
%     
%     saveas(FigS1BC, ['figs/FigS1BC.' char(ext)]);
% 
% end
