% Gillespie algorithm for Lotka Volterra dynamics
function stochastic_LV
%reactions are
% U-> 2U rate alpha
% U + V -> 2V rate beta
% V -> 0 rate gamma

global alp bet gam N vol I J IJ

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

 %reactions are
% U-> 2U rate beta
% V -> 0 rate gamma
% U + V -> 2V rate alpha

%parameters
alp = 1;
bet = 1;
gam = 0.5;
vol = 100;
 
%nullcline

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

%initial data
u0_list = [1.5,3,4];
v0_list = [1,1,1];
for j = 1:3
    
s0 = [u0_list(j),v0_list(j)]; % initial data for integration
    
[T,S] = ode23s(@deRHS,tspan,s0); 
  figure(1)
  plot( S(:,1) , S(:,2),'linewidth',2)
  xlabel(' u ','fontsize',20)
  ylabel(' v ','fontsize',20)
   
  hold on
end

figure(1)
plot([gam/alp,gam/alp],[0,5],'--', [0,5],[bet/alp,bet/alp], '--', 'linewidth',2)
 hold off

 figure(2)
 plot(T,S(:,1),T,S(:,2),'linewidth',2)
legend('U','V','fontsize',20)
 xlabel('t','fontsize',20)
  
 %now  simulate the stochastic LV birth-death process
 % set the parameters
alpha = alp; 
beta = bet/vol;
gamma = gam;

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

%Specify the change matrix for the reactions

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

%warning this can take a while if K is large
K = 2000;  % number of trials to run
 
 
%initialize the state space
x  = 10*ones(K,1);   
y = 30*ones(K,1);
 X=[];
 Y=[];
 
T = zeros(K,1);
j = 1;
ndx=find(x.*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 ;
    h(:,2) = c(2)*x.*y;
    h(:,3) = c(3)*y;
     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
     
    end
    % save one trajectory
    X(j) = x(1);
    Y(j) = y(1);
    %update ndx
    ndx=find(x.*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/vol,T(1,:),Y/vol,'linewidth',2)
 
   legend('U','V','fontsize',20)
 xlabel('t','fontsize',20)
  hold on


 hold off

hold off
figure(3)
plot(X/vol,Y/vol,'linewidth',2)
hold on
plot([gam/alp,gam/alp],[0,5],'--', [0,5],[bet/alp,bet/alp], '--', 'linewidth',2)
hold off
 xlabel(' u ','fontsize',20)
  ylabel(' v ','fontsize',20)
 
% save the time to extinction
Tex = T(:,end);
 
figure(5)

hist(Tex,50)


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

s_prime = [dudt dvdt ]';
 


