% Gillespie algorithm for a michaelis menten
% initialize the states
clear

%set the rate constants
c(1) = 2;
c(2) = 100;
c(3) = 1;

%Specify number of molecules
N = 100;
E = 50;  %enzyme molecules
 Z = zeros(N+1,E+1);
e = E;
s = N;
x(1) = s;
x(2) = e;
x(3) = 0;
C = [ -1,-1,0;1,1,0;  0,1,1];

%initialize the state space
 
X(1,:) = x; 
T(1) = 0;
j = 1;
H = 1;
K = 200000; % number of time steps
while (H>0&j<K)
    j = j+1;
    h(1) = c(1)*x(1) *x(2);
    h(2) = c(2)*(E-x(2));
    h(3) = c(3)*(E-x(2)) ;
    hc = cumsum(h);
    H = sum(h);
    % if H = 0, the population went extinct
    rn = rand(1);
    if (H>0)
    T(j) = -log(rn)/H +T(j-1);

    rn = rand(1);
    rk = min(find(rn<=hc/H)); % this determines which reaction occurs
    x = x + C(rk,:);
    X(j,:) = x;
    Z(x(1)+1,x(2)+1)=Z(x(1)+1,x(2)+1)+1;
    end
end
figure(1)
plot(X(:,1),X(:,2))
xlabel('substrate','fontsize',18)
ylabel('enzyme','fontsize',18)
 
figure(2)
mesh(Z)
figure(3)
plot(T,X(:,1)+E-X(:,2),T,X(:,1))
xlabel('time','fontsize',18)
ylabel('substrate','fontsize',18)
 