function out=gate(V,Conc,kf,kr)
%
% Compute the coefficient matrix for the system of ODEs at the potential V
% and drug concentration, Conc.  
%
% Units are /msec
%

ah = 0.126 * exp(-(V + 77.)/4.0); 
bh = 1.7 / (exp(-(V + 22.5)/12.) + 1.0);
aj = 0.055 * exp(-(V + 78.)/4.0)/(exp(-(V + 78.)/5.0) + 1.0);
bj = 0.3 / (exp(-(V + 32.)/10.) + 1.0);

k21 = ah;
k12 = bh;
k23 = Conc*4e2;			% Association rate = 4 x 10**5 /M/s
k32 = .001	;		% Dissociation rate = 1  /s
k34 = ah;			% Force inactivation to obey mass balance
k43 = bh;
k14 = Conc*kf;			% rest blocking rate
k41 = kr;			% rest unblocking rate

% uncomment the rate constant you want to force detailed balance

if Conc > 0
  %k43 = k12*k23*k34*k41 / (k21*k32*k14);		%modify beta(h) 
   k34 = k21*k32*k43*k14 / (k12*k23*k41);		%modify alpha(h) 
end

out = [[-(k12+k14) k21 0 k41]' [k12 -(k21+k23) k32 0]' ...
     [0 k23 -(k32+k34) k43]' [k14 0 k34 -(k43+k41)]']' ;