function Q = getModelQMatrix(p)

% Uses binary elements representation
% (see https://doi.org/10.1085/jgp.201411183)

% model rate constants ([/X]/s)
% ko = pore opening
% kc = pore closing
% konH = H binding rate
% koffH = H unbinding rate
% konB = Zn binding rate to blocking site
% koffB = Zn unbinding rate from blocking site
% konP = Zn binding rate to potentiating site
% koffP = Zn unbinding rate from potentiating site

ko = p(1);
kc = p(2);
konH = p(3);
koffH = p(4);
konB = p(5);
koffB = p(6);
konP = p(7);
koffP = p(8);

% Interaction energies between the pore gate (G) and the H site (H)
% H1G1 = Interaction energy for open and proton bound states (kcal/mol)
% H1Gt = Interaction energy for pore gate transition state in proton bound states (kcal/mol)
% HtG1 = Interaction energy for proton site transition state in open states (kcal/mol)

H1G1 = p(9);
H1Gt = p(10);
HtG1 = p(11);

% Interaction energies between the pore gate (G) and the Zn potentiating site (P)
% P1G1 = Interaction energy for open and potentiated states (kcal/mol)
% P1Gt = Interaction energy for pore gate transition state in potentiated states (kcal/mol)
% PtG1 = Interaction energy for potentiated site transition state in open states (kcal/mol)

P1G1 = p(12);
P1Gt = p(13);
PtG1 = p(14);

% Interaction energies between the H site (P) and the Zn potentiating site (P)
% H1P1 = Interaction energy for proton bound and potentiated states (kcal/mol)
% H1Pt = Interaction energy for potentiated site transition state in proton bound states (kcal/mol)
% HtP1 = Interaction energy for proton site transition state in potentiated states (kcal/mol)

H1P1 = p(15);
H1Pt = p(16);
HtP1 = p(17);

% pH -> H+ concentration

pH = p(18);
H = 10^(-pH); % M

% Zn concentration 

Zn = p(19); % M

% state indices are binary rep of model elements (gate,proton,blocked,potentiated)
% e.g. 0010 -> gate closed, proton site unocuppied, blocking site occupied, potentiated site unoccupied

RT = 0.593; % kcal/mol at room temp

% G H B P <-- gate, proton site, blocking site, potentiating site
Q = zeros(2^4, 2^4);

Q(1+bin2dec('0000'),1+bin2dec('1000')) = ko;
Q(1+bin2dec('0100'),1+bin2dec('1100')) = ko * exp(-H1Gt / RT);
Q(1+bin2dec('0010'),1+bin2dec('1010')) = ko;
Q(1+bin2dec('0001'),1+bin2dec('1001')) = ko * exp(-P1Gt / RT);
Q(1+bin2dec('0110'),1+bin2dec('1110')) = ko * exp(-H1Gt / RT);
Q(1+bin2dec('0101'),1+bin2dec('1101')) = ko * exp(-H1Gt / RT) * exp(-P1Gt / RT);
Q(1+bin2dec('0011'),1+bin2dec('1011')) = ko * exp(-P1Gt / RT);
Q(1+bin2dec('0111'),1+bin2dec('1111')) = ko * exp(-H1Gt / RT) * exp(-P1Gt / RT);

Q(1+bin2dec('1000'),1+bin2dec('0000')) = kc;
Q(1+bin2dec('1100'),1+bin2dec('0100')) = kc * exp(-H1Gt / RT) / exp(-H1G1 / RT);
Q(1+bin2dec('1010'),1+bin2dec('0010')) = kc;
Q(1+bin2dec('1001'),1+bin2dec('0001')) = kc * exp(-P1Gt / RT) / exp(-P1G1 / RT);
Q(1+bin2dec('1110'),1+bin2dec('0110')) = kc * exp(-H1Gt / RT) / exp(-H1G1 / RT);
Q(1+bin2dec('1101'),1+bin2dec('0101')) = kc * exp(-H1Gt / RT) / exp(-H1G1 / RT) * exp(-P1Gt / RT) / exp(-P1G1 / RT);
Q(1+bin2dec('1011'),1+bin2dec('0011')) = kc * exp(-P1Gt / RT) / exp(-P1G1 / RT);
Q(1+bin2dec('1111'),1+bin2dec('0111')) = kc * exp(-H1Gt / RT) / exp(-H1G1 / RT) * exp(-P1Gt / RT) / exp(-P1G1 / RT);

Q(1+bin2dec('0000'),1+bin2dec('0100')) = H * konH;
Q(1+bin2dec('1000'),1+bin2dec('1100')) = H * konH * exp(-HtG1 / RT);
Q(1+bin2dec('0010'),1+bin2dec('0110')) = H * konH;
Q(1+bin2dec('0001'),1+bin2dec('0101')) = H * konH * exp(-HtP1 / RT);
Q(1+bin2dec('1010'),1+bin2dec('1110')) = H * konH * exp(-HtG1 / RT);
Q(1+bin2dec('1001'),1+bin2dec('1101')) = H * konH * exp(-HtG1 / RT) * exp(-HtP1 / RT);
Q(1+bin2dec('0011'),1+bin2dec('0111')) = H * konH * exp(-HtP1 / RT);
Q(1+bin2dec('1011'),1+bin2dec('1111')) = H * konH * exp(-HtG1 / RT) * exp(-HtP1 / RT);

Q(1+bin2dec('0100'),1+bin2dec('0000')) = koffH;
Q(1+bin2dec('1100'),1+bin2dec('1000')) = koffH * exp(-HtG1 / RT) / exp(-H1G1 / RT);
Q(1+bin2dec('0110'),1+bin2dec('0010')) = koffH;
Q(1+bin2dec('0101'),1+bin2dec('0001')) = koffH * exp(-HtP1 / RT) / exp(-H1P1 / RT);
Q(1+bin2dec('1110'),1+bin2dec('1010')) = koffH * exp(-HtG1 / RT) / exp(-H1G1 / RT);
Q(1+bin2dec('1101'),1+bin2dec('1001')) = koffH * exp(-HtG1 / RT) / exp(-H1G1 / RT) * exp(-HtP1 / RT) / exp(-H1P1 / RT);
Q(1+bin2dec('0111'),1+bin2dec('0011')) = koffH * exp(-HtP1 / RT) / exp(-H1P1 / RT);
Q(1+bin2dec('1111'),1+bin2dec('1011')) = koffH * exp(-HtG1 / RT) / exp(-H1G1 / RT) * exp(-HtP1 / RT) / exp(-H1P1 / RT);

Q(1+bin2dec('0000'),1+bin2dec('0010')) = Zn * konB;
Q(1+bin2dec('1000'),1+bin2dec('1010')) = Zn * konB;
Q(1+bin2dec('0100'),1+bin2dec('0110')) = Zn * konB;
Q(1+bin2dec('0001'),1+bin2dec('0011')) = Zn * konB;
Q(1+bin2dec('1100'),1+bin2dec('1110')) = Zn * konB;
Q(1+bin2dec('1001'),1+bin2dec('1011')) = Zn * konB;
Q(1+bin2dec('0101'),1+bin2dec('0111')) = Zn * konB;
Q(1+bin2dec('1101'),1+bin2dec('1111')) = Zn * konB;

Q(1+bin2dec('0010'),1+bin2dec('0000')) = koffB;
Q(1+bin2dec('1010'),1+bin2dec('1000')) = koffB;
Q(1+bin2dec('0110'),1+bin2dec('0100')) = koffB;
Q(1+bin2dec('0011'),1+bin2dec('0001')) = koffB;
Q(1+bin2dec('1110'),1+bin2dec('1100')) = koffB;
Q(1+bin2dec('1011'),1+bin2dec('1001')) = koffB;
Q(1+bin2dec('0111'),1+bin2dec('0101')) = koffB;
Q(1+bin2dec('1111'),1+bin2dec('1101')) = koffB;

Q(1+bin2dec('0000'),1+bin2dec('0001')) = Zn * konP;
Q(1+bin2dec('1000'),1+bin2dec('1001')) = Zn * konP * exp(-PtG1 / RT);
Q(1+bin2dec('0100'),1+bin2dec('0101')) = Zn * konP * exp(-H1Pt / RT);
Q(1+bin2dec('0010'),1+bin2dec('0011')) = Zn * konP;
Q(1+bin2dec('1100'),1+bin2dec('1101')) = Zn * konP * exp(-PtG1 / RT) * exp(-H1Pt / RT);
Q(1+bin2dec('1010'),1+bin2dec('1011')) = Zn * konP * exp(-PtG1 / RT);
Q(1+bin2dec('0110'),1+bin2dec('0111')) = Zn * konP * exp(-H1Pt / RT);
Q(1+bin2dec('1110'),1+bin2dec('1111')) = Zn * konP * exp(-PtG1 / RT) * exp(-H1Pt / RT);

Q(1+bin2dec('0001'),1+bin2dec('0000')) = koffP;
Q(1+bin2dec('1001'),1+bin2dec('1000')) = koffP * exp(-PtG1 / RT) / exp(-P1G1 / RT);
Q(1+bin2dec('0101'),1+bin2dec('0100')) = koffP * exp(-H1Pt / RT) / exp(-H1P1 / RT);
Q(1+bin2dec('0011'),1+bin2dec('0010')) = koffP;
Q(1+bin2dec('1101'),1+bin2dec('1100')) = koffP * exp(-PtG1 / RT) / exp(-P1G1 / RT) * exp(-H1Pt / RT) / exp(-H1P1 / RT);
Q(1+bin2dec('1011'),1+bin2dec('1010')) = koffP * exp(-PtG1 / RT) / exp(-P1G1 / RT);
Q(1+bin2dec('0111'),1+bin2dec('0110')) = koffP * exp(-H1Pt / RT) / exp(-H1P1 / RT);
Q(1+bin2dec('1111'),1+bin2dec('1110')) = koffP * exp(-PtG1 / RT) / exp(-P1G1 / RT) * exp(-H1Pt / RT) / exp(-H1P1 / RT);

% unitary matrix (neg sum or each row on the diagonal)
Q = Q - diag(sum(Q,2));

end