% Gillespie simulation for a birth-death process with bimolecular death
% reaction 1 death S +S -> S+I rate alpha
% reaction 2 birth S -> 2S  rate beta
function birth_death_stochastic
% initialize the states
global Nde  alpha beta
 N = 10; %initial population size
Nde = 30;  %this is the number of p_k's to simulate, needed because of birth
K = 10000;  % number of trials
alpha = 1;  % without loss of generality
beta = 5;  %  

%set the rate constants for the two reactions:
c(1) = alpha;  %S +S -> S+ I  % I tracks the number of deaths
c(2) = beta; %S -> 2S

%Specify the change matrix

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

%initialize the state space

s = N* ones(K,1);
i = zeros(K,1); %start each population with no deaths yet
  
S = s; %This will keep track of the full history/trajectory
I = i;
T = zeros(K,1);  %This will track the transition times
j = 1;
rk = [];
%the population cannot go extinct

Kt = 500; % max number of time steps
while (j<Kt)
    j = j+1;
    h(:,1) = c(1)*s.*(s-1)/2;  %reaction rate 1
    h(:,2) = c(2)*s; % reaction rate 2
    hc = cumsum(h')';
    H = sum(h')';
    % if H = 0, i is zero no reaction is possible: we're done
    rn = rand(K,2);
    
    Nt = sum(H>0);
    T(:,j) = T(:,j-1);  % add the current time to T
    for k = 1:K
    if (H(k)>0)
    T(k,j) = - log(rn(k,1))/H(k) +T(k,j-1); % time of next reaction
    rk  = min(find(rn(k,2) <=hc(k,:)/H(k) )); % this determines which reaction occurs
    s(k) = s(k) + C(rk,1); % update s 
    i(k) = i(k) + C(rk,2);  % update i
    
    end
    S(:,j) = s; % store the current state
    I(:,j) = i;
    end
end
 
nj=j;  %The number of times step taken
 kj = 1;; % use the first trial to plot a sample trajectory 
 tt = reshape([T(kj,:);T(kj,:)],2*nj,1);  % make a plot with stepwise increments
 nn = reshape([S(kj,:);S(kj,:)],2*nj,1); 
  
%plot one trajectory
 figure(1)
 plot(tt(2:end),nn(1:end-1,1),'linewidth',2)
 axis([0 5 0 Nde])
 xlabel('t', 'fontsize', 20)
 ylabel('Particle Number', 'fontsize', 20)
 
%
 % find the distribution at some fixed time
 Tm = min(T(:,end));
 neq= zeros(Nde,1);
 
 for j = 1:K
     indx = min(find(T(j,:)>=Tm)); % this is the index with time Tm
     eq = S(j,indx);
     
     % the state at this time
     neq(eq) = neq(eq) + 1; %create a histogram of the values of S
 end
 
% simulate the de's
p0 = zeros(Nde,1); %initialize p_k 
p0(N) = 1;
 
tstep = 0.1;
t_end = 2;  %check to make sure p_k has reached equilibrium
tspan = [0:tstep:t_end];   
[T,S] = ode23s(@deRHS,tspan,p0); 
 
figure(2)
for j = 1:Nde
plot(T,S(:,j),'linewidth',2)
hold on
end
hold off
xlabel('t','fontsize',20)
ylabel('p_k(t)','fontsize',20)
 
figure(3)
plot([1:Nde],neq/K,'*',[1:Nde],S(end,:),'--','linewidth',2)
xlabel('n','fontsize',20)
ylabel('% at level n','fontsize',20)

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


