function [time, I, Po, PG, PH, PB, PP, idx] = supfig1_simulation(p, dt, Zn, washoutSec)

% p = [ko, kc, konB, koffB, konP, koffP, P1G1, P1Gt, PtG1]
% dt (sec)
% Zn (M)

% ko = pore opening (/M/s)
% kc = pore closing (/s)
% konB = Zn binding rate to blocking site (/M/s)
% koffB = Zn unbinding rate from blocking site (/s)
% konP = Zn binding rate to potentiating site (/M/s)
% koffP = Zn unbinding rate from potentiating site (/s)
% P1G1 = Interaction energy stabilizing open pore gate and Zn potentiated states (kcal/mol)
% P1Gt = Interaction energy for pore gate transition state in Zn potentiated states (kcal/mol)
% PtG1 = Interaction energy for Zn potentiation site transition state in open pore gate states (kcal/mol)

% pH 7.4 -> pH 7.4 + Zn -> pH 5.5 -> pH 7 .4

wash = washoutSec;%[0, 1, 4, 8, 16]; % Zn washout duration (s)

time = (0:dt:2+max(wash)+16+8+4)'; % s

N = 2^4; % # of states
K = length(wash); % # sweeps

P = zeros(N,length(time),K); % state probability arrays (state,t,sweep)

for i = 1:K
    clear Q
    if wash(i) == 0
        t{1} = 0:dt:2+max(wash); % pH 7.4
        t{2} = 2+max(wash):dt:2+max(wash)+16; % pH 7.4 + Zn
        t{3} = 2+max(wash)+16:dt:2+max(wash)+16+8; % pH 5.5
        t{4} = 2+max(wash)+16+8:dt:2+max(wash)+16+8+4; % pH 7.4

        idx{1} = find(time <= t{1}(end));
        idx{2} = find((time >= t{2}(1)) & (time <= t{2}(end)));
        idx{3} = find((time >= t{3}(1)) & (time <= t{3}(end)));
        idx{4} = find(time >= t{4}(1));

        Q{1} = getModelQMatrix([p,7.4,0]);
        Q{2} = getModelQMatrix([p,7.4,Zn]);
        Q{3} = getModelQMatrix([p,5.5,0]);
        Q{4} = getModelQMatrix([p,7.4,0]);
    else
        t{1} = 0:dt:2+max(wash)-wash(i); % pH 7.4
        t{2} = 2+max(wash)-wash(i):dt:2+max(wash)-wash(i)+16; % pH 7.4 + Zn
        t{3} = 2+max(wash)-wash(i)+16:dt:2+max(wash)+16; % pH 7.4
        t{4} = 2+max(wash)+16:dt:2+max(wash)+16+8; % pH 5.5
        t{5} = 2+max(wash)+16+8:dt:2+max(wash)+16+8+4; % pH 7.4

        idx{1} = find(time <= t{1}(end));
        idx{2} = find((time >= t{2}(1)) & (time <= t{2}(end)));
        idx{3} = find((time >= t{3}(1)) & (time <= t{3}(end)));
        idx{4} = find((time >= t{4}(1)) & (time <= t{4}(end)));
        idx{5} = find(time >= t{5}(1));

        Q{1} = getModelQMatrix([p,7.4,0]);
        Q{2} = getModelQMatrix([p,7.4,Zn]);
        Q{3} = getModelQMatrix([p,7.4,0]);
        Q{4} = getModelQMatrix([p,5.5,0]);
        Q{5} = getModelQMatrix([p,7.4,0]);
    end
    
    % first epoch state probabilites at equilibrium
    S = [Q{1} ones(N,1)];
    u = ones(1,N);
    P0 = u * inv(S * (S'));
    P0(P0 < 0) = 0;
    P0 = P0 ./ sum(P0);
    P(:,idx{1},i) = repmat(P0',1,length(idx{1}));
    
    % remaining epoch state probabilities
    for j = 2:length(Q)
        timej = t{j} - t{j}(1);
        for k = 2:length(t{j})
            P(:,idx{j}(k),i) = P0 * expm(Q{j} .* timej(k));
        end
        P0 = P0 * expm(Q{j} .* timej(end));
    end
end

Ostates = [1+bin2dec('1000'), 1+bin2dec('1100'), 1+bin2dec('1001'), 1+bin2dec('1101')];
Gstates = [1+bin2dec('1000'), 1+bin2dec('1100'), 1+bin2dec('1010'), 1+bin2dec('1001'), 1+bin2dec('1110'), 1+bin2dec('1101'), 1+bin2dec('1011'), 1+bin2dec('1111')];
Hstates = [1+bin2dec('0100'), 1+bin2dec('1100'), 1+bin2dec('0110'), 1+bin2dec('0101'), 1+bin2dec('1110'), 1+bin2dec('1101'), 1+bin2dec('0111'), 1+bin2dec('1111')];
Bstates = [1+bin2dec('0010'), 1+bin2dec('1010'), 1+bin2dec('0110'), 1+bin2dec('0011'), 1+bin2dec('1110'), 1+bin2dec('1011'), 1+bin2dec('0111'), 1+bin2dec('1111')];
Pstates = [1+bin2dec('0001'), 1+bin2dec('1001'), 1+bin2dec('0101'), 1+bin2dec('0011'), 1+bin2dec('1101'), 1+bin2dec('1011'), 1+bin2dec('0111'), 1+bin2dec('1111')];

Po = squeeze(sum(P(Ostates,:,:), 1));
PG = squeeze(sum(P(Gstates,:,:), 1));
PH = squeeze(sum(P(Hstates,:,:), 1));
PB = squeeze(sum(P(Bstates,:,:), 1));
PP = squeeze(sum(P(Pstates,:,:), 1));

% ??? Nernst reversal does not seem to match data (see Fig. 5 in first paper)
% EH74 = 60 * log10(10^-7.4 / 10^-7.5); % mV
% EH55 = 60 * log10(10^-5.5 / 10^-7.5); % mV

% ??? approx reversal as -80 mV at pH 7.4 and 0 mV at pH 5.5
EH = repmat(-80, size(time)); % mV
EH((time > 2+max(wash)+16) & (time <= 2+max(wash)+16+8)) = 0;
Vhold = -80; % mV
conductance = 1; % pS

I = Po .* repmat((Vhold - EH) * 1e-3 * conductance, [1,size(Po,2)]); % pA

end