% Gillespie simulation for the enzyme decay process
% reaction 1 complex formation S + E -> C  rate alpha
% reaction 2 complex degradation C -> S + E  rate beta
% product formation C-> P + E  rate gamma

function enzymatic_decay
global N alpha beta gamma
N = 40; %initial substrate population size
% enzyme is initially 1;
% complex is initially zero
%
K = 2500;  % number of trials
alpha = 1;  % without loss of generality
beta = 2;   
gamma = 1;
%set the rate constants for three reactions:
c(1) = alpha; %S + E -> C  rate alpha%  
c(2) = beta; %C -> S + E  rate beta
c(3) = gamma; %C-> P + E  rate gamma
%Specify the change matrix

C = [-1,-1,1;1,1,-1;0,1,-1];  %What happens to s, e and c when a reaction occurs:

%initialize the state space

s = N*ones(K,1);  %N substrate molecules
e = ones(K,1); %start each population with one enzyme molecule
cm = zeros(K,1);  % start with zero complex molecules
S = s; %This will keep track of the full history/trajectory
E = e;
Cm = cm;
T = zeros(K,1);  %This will track the transition times
j = 1; % j keeps track of the number of time steps
rk = [];
Nt = K;  % number of trials that have not yet gone extinct; extinction means s=0
ndx = find(s>0);

Kt = 5000; % max number of time steps
while (Nt>0&j<Kt)
    j = j+1;
    h(:,1) = c(1)*s.*e;  %reaction rate 1
    h(:,2) = c(2)*cm; % reaction rate 2
    h(:,3) = c(3)*cm; % reaction rate 3
    hc = cumsum(h')';
    H = sum(h')';
    % if H = 0,  no reaction is possible: we're done
     
    
     rn = rand(Nt,2); %find random number for the number of not yet terminated trajectories
    
    T(:,j) = T(:,j-1);  % add the current time to T
    
    for k = 1:Nt
     % update only those for which H >0
     
    T(ndx(k),j) = - log(rn(k,1))/H(ndx(k)) +T(ndx(k),j-1); % time of next reaction
    rk  = min(find(rn(k,2) <=hc(ndx(k),:)/H(ndx(k)) )); % this determines which reaction occurs
    
    s(ndx(k)) = s(ndx(k)) + C(rk,1); % update s 
    e(ndx(k)) = e(ndx(k)) + C(rk,2); % update e
    cm(ndx(k)) = cm(ndx(k)) + C(rk,3); % update cm
    
    end
    S(:,j) = s;  % store the current states
    E(:,j) = e;
    Cm(:,j) = cm;
    ndx = find(s+cm>0);  %check to see which trajectories are not extinct yet
    if (isempty(ndx)==1)
        Nt = 0;
    else
        Nt = length(ndx); %Nt is the number of trajectories not yet extinct
    end
end
nj = j; % maximum number of reactions

 sj= 1; % use the first  trial to plot a sample trajectory 
tt = reshape([T(sj,:);T(sj,:)],2*nj,1);  % make a plot with stepwise decay
ss = reshape([S(sj,:);S(sj,:)],2*nj,1);
%plot one decay trajectory
 figure(1)
 plot(tt(1:end-1),ss(2:end),'linewidth',2)
 xlabel('t', 'fontsize', 20)
 ylabel('Substrate Number', 'fontsize', 20)  
 
  
% create the pdf for extinction times from the data
 
 
 [NN,TT]=hist(T(:,end),50);  % histogram of the extinction times with 50 boxes
 tt = sort(TT);
   
 dt = mean(TT(2:end)-TT(1:end-1));  %timestep increment on the histogram
 
 % NN/(K*dt) is the approximate pdf for the extinction times from data
 % it is normalized to have total integral = 1
  Mt = mean(T(:,end))
  Vt = var(T(:,end))
  figure(2)
 plot(TT, NN/(K*dt),'*',TT,exp(-(TT-Mt).^2/(2*Vt))/(sqrt(2*pi*Vt)),'--','linewidth',2)
 xlabel('extinction time','fontsize',20)
 ylabel('pdf of extinction times','fontsize',20)

  
  % numerical solution of the master equations
  % initialize the probabilities
p0 = [zeros(2*N+1,1);1]; %initialize p_{k,j} 

 tstep = 0.02;
t_end = 100;
tspan = [0:tstep:t_end];   
[T,S] = ode23s(@deRHS,tspan,p0); 
 figure(3)
 plot(T,S(:,2))
   
 figure(2)
 plot(TT, NN/(K*dt),'*',TT,exp(-(TT-Mt).^2/(2*Vt))/(sqrt(2*pi*Vt)),'r--', ...
     T,gamma* S(:,1),'k--','linewidth',2)
 xlabel('extinction time','fontsize',20)
 ylabel('pdf of extinction times','fontsize',20)
 axis([0 100 0 0.05])
function p_prime=deRHS(t,p)  % right hand side for ode system
global N alpha beta gamma
  
K = [0:N]';
p0 = p([1:N+1]);
p1 = p([N+2:2*N+2]);

p1_prime = -(alpha*K).*p1+beta*[0;p0(1:end-1)]+gamma*p0;
p0_prime = -(beta+gamma)*p0 + alpha*[K(2:end);0].*[p1(2:end);0];

p_prime = [p0_prime;p1_prime];


         