function Q = Q(V, concentration, drug)
% P' = [O C1 C2 C3 I I1 I2 I3 ID I1D I2D I3D]

% Voltage rates
[am bm ah bh] = HHrates(V);

% Drug rates
switch drug
	case 'phenytoin'
		bd = 10*concentration;
		ad = 7e-6*10;
	case 'carbamazepine'
		bd = 38*concentration;
		ad = 25e-6*38;
	otherwise
		bd = 0;
		ad = 1; % guarantee unique stable state
end

% Q(i,j) = P(i)->P(j)
%   O   C1   C2  C3   I   I1   I2  I3  ID  I1D  I2D  I3D
Q = [ ...
	0 3*bm    0   0  bh    0    0   0   0    0    0    0; ...
   am    0 2*bm   0   0   bh    0   0   0    0    0    0; ...
	0 2*am    0  bm   0    0   bh   0   0    0    0    0; ...
	0    0 3*am   0   0    0    0  bh   0    0    0    0; ...
   ah    0    0   0   0 3*bm    0   0  bd    0    0    0; ...
	0   ah    0   0  am    0 2*bm   0   0   bd    0    0; ...
	0    0   ah   0   0 2*am    0  bm   0    0   bd    0; ...
	0    0    0  ah   0    0 3*am   0   0    0    0   bd; ...
	0    0    0   0  ad    0    0   0   0 3*bm    0    0; ...
	0    0    0   0   0   ad    0   0  am    0 2*bm    0; ...
	0    0    0   0   0    0   ad   0   0 2*am    0   bm; ...
	0    0    0   0   0    0    0  ad   0    0 3*am    0; ...
	];

%    O   C1   C2  C3   I   I1   I2  I3  ID
% Q = [ ...
% 	0 3*bm    0   0  bh    0    0   0   0; ...
% 	am    0 2*bm   0   0   bh    0   0   0; ...
% 	0 2*am    0  bm   0    0   bh   0   0; ...
% 	0    0 3*am   0   0    0    0  bh   0; ...
% 	ah    0    0   0   0 3*bm    0   0  bd; ...
% 	0   ah    0   0  am    0 2*bm   0   0; ...
% 	0    0   ah   0   0 2*am    0  bm   0; ...
% 	0    0    0  ah   0    0 3*am   0   0; ...
% 	0    0    0   0  ad    0    0   0   0; ...
% 	];

Q = Q - diag(sum(Q,2));