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

% 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 5.5 -> pH 5.5 + Zn -> pH 5.5 -> pH 7 .4

t{1} = 0:dt:1; % pH 7.4
t{2} = 1:dt:5.1; % pH 5.5
t{3} = 5.1:dt:7; % pH 5.5 + Zn
t{4} = 7:dt:11; % pH 5.5
t{5} = 11:dt:12; % pH 7.4

time = (0:dt:12)';
idx{1} = find(time <= t{1}(end)+dt/2);
idx{2} = find((time >= t{2}(1)-dt/2) & (time <= t{2}(end)+dt/2));
idx{3} = find((time >= t{3}(1)-dt/2) & (time <= t{3}(end)+dt/2));
idx{4} = find((time >= t{4}(1)-dt/2) & (time <= t{4}(end)+dt/2));
idx{5} = find(time >= t{5}(1)-dt/2);

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

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

for i = 1:K
    Q{1} = getModelQMatrix([p,7.4,0]);
    Q{2} = getModelQMatrix([p,5.5,0]);
    Q{3} = getModelQMatrix([p,5.5,Zn(i)]);
    Q{4} = getModelQMatrix([p,5.5,0]);
    Q{5} = getModelQMatrix([p,7.4,0]);
    
    % 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:5
        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
for j = 2:4
    EH(idx{j}) = 0;
end
Vhold = -80; % mV
conductance = 1; % pS

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


end