function    Chem_React 
% This code simulates the Chemical Reaction on 
% page 2350 (Eqs.38a-38c) of Gillespie's Paper 
% using the Stochastic Queuing routine from Figure 2 
% of the same paper. We recreate the images in 
% Figures 8a and 8c (not exactly the same since 
% simulations are intrinsicly stochastic).
%
%       Written 11-24-2007 by John L. Aven

clc
clear all
close all

% -----      Initialize Molecular Populations  
Y_1_i = 1000.0;
Y_2_i = 1000.0;

% -----   FoodStuff Population(Quantity)  
X = 10.0 ^ (1.0);

% -----   The Number of Reactions  
M = 3;

% -----  Initialize Chemical Reaction Coefficients 
c1 = 10.0 / X;
c2 = 0.01;
c3 = 10.0;

%   --   Place coefficients in Vector C 
C = [c1;c2;c3];

tau = 0.1;                  %The Waiting Time.

% Initialization of Evolutionary Variables;
Y = [Y_1_i;Y_1_i];
time = 0.0;
n = 1;
t = 0;
t_max = 10.0;

while time <= t_max   
%     C(1) = 10.0/X;     % Uncomment if you allow X
                         % to change with time.
    A = C .* Reaction_Combination_Functions([X;Y(:,n)]);
    
    A_sum = sum(A);         % COEFFICIENT a_o
    
    u = rand(2);            % The two uniform (0,1) random variables.
    
    tau = ( 1 / A_sum ) * log ( 1 / u(1) );     % The next waiting time.
                    
    time = time + tau;                          % Update the time
    
    t(n+1) = time ;                             % Store time for plotting
    
    seek_value = u(2) * A_sum;                  % next reaction threshold
    
    %  ---  Find the Reaction to Enact  --  %
    seek_truth = 1;
    seek_test  = A(1);
    seek_mu    = 1;
    while (seek_truth == 1)
        
        if (seek_value < seek_test)           
            mu = seek_mu;
            seek_truth = 0;
            
        elseif (seek_mu < M)          
            seek_mu = seek_mu + 1;
            seek_test = seek_test + A(seek_mu);
            
        else        
            mu = M;
            seek_truth = 0;
            
        end        
    end
    
    n = n + 1;                              % Update the reaction counter    
    if (mu == 1)        % Reaction 1 (X + Y1 --> 2*Y1)      
        Y(1,n) = Y(1,n-1) + 1.0;        
%         X = X - 1.0;  % Uncomment if X can change
                        % with time
        Y(2,n) = Y(2,n-1);
        
    elseif (mu == 2)    % Reaction 2 (Y1 + Y2 --> 2*Y2)       
        Y(1,n) = Y(1,n-1) - 1.0;
        Y(2,n) = Y(2,n-1) + 1.0;
        
    elseif (mu == 3)    % Reaction 2 (Y2 --> Z)       
        Y(1,n) = Y(1,n-1);
        Y(2,n) = Y(2,n-1) - 1.0; 
    end   
end

figure(1)

plot(t,Y(1,:),'k')
axis tight
xlabel('Time','FontSize',20)
ylabel('Y_1','FontSize',20)

figure(2)

plot(t,Y(2,:),'k')
axis tight
xlabel('Time','FontSize',20)
ylabel('Y_2','FontSize',20)

figure(3)

plot(t,Y(1,:),'k')
hold on
plot(t,Y(2,:),'r')
axis tight
xlabel('Time','FontSize',20)
ylabel('Y_i','FontSize',20)
legend('Y_1','Y_2')

figure(4)

plot(Y(1,:),Y(2,:),'k.','MarkerSize',1)
axis tight
xlabel('Y_1','FontSize',20)
ylabel('Y_2','FontSize',20)
title('Dot Plot')

figure(5)

plot(Y(1,:),Y(2,:),'k','MarkerSize',1)
axis tight
xlabel('Y_1','FontSize',20)
ylabel('Y_2','FontSize',20)
title('Continuous Plot')

size(Y)

% ----------------------------------------
function [H] = Reaction_Combination_Functions(Var)
% The Reaction Function for all Reactions in a Vectorized Form.

Var2 = [Var(2:3);1.0];
H = Var .* Var2;