% Gillespie algorithm for SIRS dynamics
function stochastic_SIRS
%reactions are
% S+I->2I rate alpha
% I->R rate beta
% R->S rate gamma

global alp bet gam s0

 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

%parameters
alp = 1.5;
bet = 1.2;
gam = 2.1;
vol = 100;
s0=1;

 
%nullclines


figure(1)
s = [0:.01:s0];
i0 = gam*(s0-s)./(alp*s+gam);
str = bet/alp;
istr = gam*(s0-str)./(alp*str+gam)
plot(s,i0,'--',bet*ones(1,length(s))/alp,s,'--',str,istr,'k*',1,0,'k*','linewidth',2)
xlabel('s/s_0','fontsize',20)
ylabel('i/s_0','fontsize',20)
hold on

% now create some phase portraits for the deterministic system
tstep = 0.1; %time step size
t_end = 25; % length of time interval
tspan = [0:tstep:t_end];

%initial data
u0_list = [0.9];
v0_list = [0.01];
for j = 1:1
    
u0v0 = [u0_list(j),v0_list(j)]; % initial data for integration
    
[T,S] = ode23s(@deRHS,tspan,u0v0); 
  figure(1)
  plot( S(:,1) , S(:,2),'linewidth',2)
   
  annotation('textarrow',[0.45,0.45],[0.55 0.45],'linewidth',2) 
   annotation('textarrow',[0.8 0.72],[0.45,0.45],'linewidth',2) 
  hold on
end

  
 %now  simulate the stochastic LV birth-death process
 % set the parameters
alpha = alp/vol; 
beta = bet;
gamma = gam;

%set the rate constants
c(1) = alpha;
c(2) = beta;
c(3) = gamma;

%Specify the change matrix for the reactions

C = [-1,1,0;0,-1,1;1,0,-1];

%warning this can take a while if K is large
K = 1000;  % number of trials to run
 
 
%initialize the state space
x  =  ones(K,1);   
y = 99*ones(K,1);
z = zeros(K,1);
 X=[];
 Y=[];
 Z = [];
 Y(1) = 99;
 X(1) = 1;
T = zeros(K,1);
j = 1;
ndx=find(y>0);
Nt = K; 
Kt = 10000; % max number of time steps, just in case it doesn't go extinct
while (Nt>0&j<Kt)
    j = j+1;
    h(:,1) = c(1)*x.*y ;
    h(:,2) = c(2)*y;
    h(:,3) = c(3)*z;
     hc = cumsum(h')';
    H = sum(h')';
     
    rn = rand(Nt,2);
    T(:,j) = T(:,j-1);  % add the current time to T
    
    for k = 1:Nt
     
    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
    
    x(ndx(k)) = x(ndx(k)) + C(rk,1); % update x
    y(ndx(k)) = y(ndx(k)) + C(rk,2); % update y
    z(ndx(k)) = z(ndx(k)) + C(rk,3); % update z
     
    end
    % save one trajectory
    X(j) = x(1);
    Y(j) = y(1);
    Z(j) = z(1);
    %update ndx
    ndx=find( y>0);
    if (isempty(ndx)==1)
        Nt = 0;
    else
        Nt = length(ndx); %Nt is the number of trajectories not yet extinct
    end
end
%plot a sample trajectory 
figure(4)
plot(T(1,:),X ,T(1,:),Y ,T(1,:),Z ,'linewidth',2)
 
   legend('S','I','R','fontsize',20)
 xlabel('t','fontsize',20)
 
figure(3)
plot(X/vol,Y/vol,'linewidth',2)
hold on
plot(s,i0,'--',bet*ones(1,length(s))/alp,s,'--','linewidth',2)
hold off
 xlabel(' s ','fontsize',20)
  ylabel(' i ','fontsize',20)
 
% save the time to extinction
Tex = T(:,end);
 
figure(5)

hist(Tex,50)
 xlabel(' t ','fontsize',20)
  ylabel('Extinction Histogram','fontsize',20)

function s_prime=deRHS(t,s)  % right hand side for ode system
global alp bet gam s0
 
u = s(1);
i = s(2);
 
dudt = -alp*u*i+gam*(s0-u-i)
dvdt = alp*u*i-bet*i
 

s_prime = [dudt dvdt ]';
 


