% Gillespie simulation for a birth-death process (Keizer's paradox)

% reaction 1 birth S + U  -> 2U rate alpha
% reaction 2 death 2U -> U  rate beta
% reaction 3 death U -> P  rate gamma
function birth_death_stochastic

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);
% initialize the states
global Nde  alpha beta gamma % N = total number of particles


alpha = 1;  % without loss of generality
beta = 0.05;  %  
gamma = 0.5;  % must have alpha > gamma

% the stable steady state is at (alpha-gamma)/beta

N = fix((alpha-gamma)/beta) %initial population size
Nde = 4*N%  We must allow for more than the initial population size 
K = 2000;  % number of trials

%set the rate constants for three reactions:
c(1) = alpha;  %S + U  -> 2U rate alpha birth S held constant
c(2) = beta; %2U -> U  rate beta  death
c(3) = gamma;  %death U -> P  rate gamma
%Specify the change matrix

C = [1,0;-1,0;-1,1];  %What happens to u and p when a reaction occurs:

%initialize the state space

u = N* ones(K,1);
p = zeros(K,1); %start each population with no product
 
U = u;
P = p;
T = zeros(K,1);  %This will track the transition times
j = 1;
rk = [];

ndx = find(u>0);  %update only the s's that are not zero
Nt = length(ndx); %=K at first
Kt = 5000; % max number of time steps to prevent very long runs without termination
while (Nt>0&j<Kt)
    j = j+1;
    h(:,1) = c(1)*u;  %reaction rate 1
    h(:,2) = c(2)*u.*(u-1); % reaction rate 2  (a quadratic rate = bimolecular reaction)
    h(:, 3) = c(3)*u;
    hc = cumsum(h')'; % the cumulative sum of h
    H = sum(h')';
   
    rn = rand(Nt,2);
    
    T(:,j) = T(:,j-1);  % add the current time to T
    
    for k = 1:Nt
     % update only those for which s >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
    u(ndx(k)) = u(ndx(k)) + C(rk,1); % update s and i
    p(ndx(k)) = p(ndx(k)) + C(rk,2);
    
    end 
    % save the values of the   trajectories
    U(:,j) = u ;
    P(:,j) = p ;
    ndx = find(u>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
 
  
 kj = 1; % use the first several trials to plot a sample trajectory 
  
% phase portrait of some trajectories
  
 
%plot one decay trajectory
 figure(1)
 plot(T(kj,:),U(kj,:),'linewidth',2)
 xlabel('t', 'fontsize', 20)
 ylabel('Particle Number', 'fontsize', 20)
 set(gca, 'YTick', 0:4:24);
  
title('Sample Trajectory','fontsize',20)
 
 % create a histogram of product
 [NN,PP]=hist(P(:,end),50);  % histogram of the final product with 50 boxes
 figure(4)
 plot(NN,PP/K,'*','linewidth',2)
 xlabel('Product', 'fontsize', 20)
 ylabel('Product histogram', '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
 
% simulate the de's
p0 = zeros(Nde+1,1); %initialize p_k 
p0(N+1) = 1;
 
 tstep = 0.1;
t_end =  max(T(:,end));
tspan = [0:tstep:t_end];   
[T,S] = ode23s(@deRHS,tspan,p0); 
 
figure(3)
for j = 1:Nde
plot(T,S(:,j))
hold on
end
xlabel('t','fontsize',20)
ylabel('p_n(t)','fontsize',20)

hold off
 
figure(2)
 plot(TT,NN/(K*dt),'*',T,gamma*S(:,2),'--','linewidth',2)
  xlabel('t','fontsize',20)
ylabel('dp_0/dt','fontsize',20)

 
 
function p_prime=deRHS(t,p)  % right hand side for ode system
global Nde alpha beta gamma
 
K = [0:Nde]'; 
 
p_prime = -((alpha+gamma)*K+beta*K.*(K-1)).*p...
+alpha*[0;K(1:end-1)].*[0;p(1:end-1)]...
    +beta*[K(2:end);0].*([K(2:end);0]-1).*[p(2:end);0] ...
    + gamma*[K(2:end);0].*[p(2:end);0];

 