% SI decay of particles via Gillespie algorithm
% the reaction is S+I->2I  with S+I constant
function stochastic_decay
global N alpha

%set the initial number of particles
N = 25; % initial particle number
%start with N-1 S and 1 I particles
% set the death rate alpha
alpha = 1;  % rate of decay
 
K = 10000; %total number of trials
 
% start with 
t = zeros(N,K);% this is where we record the reaction times
R = rand(N-1,K); % generate the needed random numbers
for j = 1:N-1
    n = N-j;% number of S particles left
    rxn = alpha*n*(N-n); % reaction rate with n particles left
    dt = -log(R(j,:))/rxn; % next decay time interval
    t(j+1,:) = t(j,:) + dt; %update the time
    
    end
    s = 1; % use the first trial for plots
 tt = reshape([t(:,s)';t(:,s)'],2*(N),1);  % make a plot with stepwise decay
  nn = reshape([ N-1:-1:0; N-1:-1:0],2*(N),1);
%plot one decay trajectory
 figure(1)
 plot(tt(2:end),nn(1:end-1,1),tt,N*(N-1)*exp(-N*tt)./(1+(N-1)*exp(-N*tt)),'--','linewidth',2)
 xlabel('t', 'fontsize', 20)
 ylabel('Particle Number', 'fontsize', 20)
 
 %plot one decay trajectory on a log scale
figure(2)
 semilogy(tt,nn./(N-nn),tt,(N-1)*exp(-N*tt),'--','linewidth',2)
 xlabel('t', 'fontsize', 20)
 ylabel('Particle Number', 'fontsize', 20)
 title('semilog plot: Particle number vs. time','fontsize',20)
 
% create the pdf for extinction times from the data
  
 [NN,TT]=hist(t(N,:),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
 
 % the sample mean:
 mean(t(N,:))
  
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % now simulate the differential equation system
 
 %dp_k/dt = a*(k+1)*(N-k-1)*p_(k+1) -a*k*(N-k)*p_k

% initialize the probabilities
p0 = [zeros(N-1,1);1]; %initialize p_k 
 tstep = 0.02;
t_end = 0.8;
tspan = [0:tstep:t_end];   
[T,S] = ode23s(@deRHS,tspan,p0); 
 
figure(3)
for j = 1:N
plot(T,S(:,j),'linewidth',2)
hold on
end
xlabel('t','fontsize',20)
ylabel('p_k(t)','fontsize',20)
hold off

% compare the data with the extinction time distribution dp_0/dt:
figure(4)
plot(TT,NN/(K*dt),'*',T,alpha*(N-1)*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 N alpha 
 
K = [0:N-1]';
NmK = N-K;
KNmK = K.*NmK;
p_prime = -alpha*KNmK.*p+alpha*[KNmK(2:end);0].*[p(2:end);0];

 

 
 
    