%
% Chemical Equations:
%
% 1) mRNAa + rRNA -> mRNAa + rRNA + A       (k_xlate)
% 2) mRNAa        -> decay                  (k_dec_mrna_a)
% 3) A + A        -> A2                     (k_dimer_a)
% 4) A            ->                        (k_dec_a)
% 5) A2           -> A + A                  (k_sngl_a)
% 6) A2           ->                        (k_dec_a2)
% 7) Pz + A2      -> PzA2                   (k_rprs_a2)
% 8) PzA2         -> Pz + A2                (k_dissoc_a2)
% 9) PzA2         -> Pz                     (k_dec_a2)
% 10) PzA4        -> PzA2                   (k_dec_a2)
% 11) PzA2 + A2   -> PzA4                   (k_rprs_a4)
% 12) PzA4        -> PzA2 + A2              (k_dissoc_a4)
% 13) Pz   + RNAp -> Pz + RNAp + mRNAz      (k_xcribe)
% 14) mRNAz       ->                        (k_dec_mrna_z)

% 15) l1 -> l1 + mRNAa                      (k_l1)

%%%
%%% Definitions of global variables shared with the other files
%%%

% Variables
global mRNAa;
global A;
global A2;
global Pz;
global PzA2;
global PzA4;
global mRNAz;

global l1;
% global l2;

% Rate constants
global k_xlate;
global k_xcribe;
global k_dec_mrna_a;
global k_dec_mrna_z;
global k_dec_a;
global k_dec_a2;
global k_dimer_a;
global k_sngl_a;
global k_rprs_a2;
global k_rprs_a4;
global k_dissoc_a2;
global k_dissoc_a4;

global k_l1;
% global k_l2;

% Various simulation definitions
global REGRESSION;         % Monte Carlo Simulation Times
global END_TIME;
global RECORD_STEP;
global EQ_NUMBERS;
global AVALS;
global PLASMID_COPIES;
global rRNA;
global RNAp;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization & Starting The Simulation 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
REGRESSION     = 1;
END_TIME       = 32400;  % 9 hours
RECORD_STEP    = 1;
% # of equations
EQ_NUMBERS     = 15;
PLASMID_COPIES = 20;
AVALS          = zeros(EQ_NUMBERS, 1);

% arrays to store current simulation values
REC_mRNAa = zeros(END_TIME/RECORD_STEP + 1, 1);
REC_A     = zeros(END_TIME/RECORD_STEP + 1, 1);
REC_A2    = zeros(END_TIME/RECORD_STEP + 1, 1);
REC_Pz    = zeros(END_TIME/RECORD_STEP + 1, 1);
REC_PzA2  = zeros(END_TIME/RECORD_STEP + 1, 1);
REC_PzA4  = zeros(END_TIME/RECORD_STEP + 1, 1);
REC_mRNAz = zeros(END_TIME/RECORD_STEP + 1, 1);
REC_time  = zeros(END_TIME/RECORD_STEP + 1, 1);

REC_l1  = zeros(END_TIME/RECORD_STEP + 1, 1);

% array to store statistics about multiple simulations
sum_REC_mRNAa = zeros(END_TIME/RECORD_STEP + 1, 1);
sum_REC_A     = zeros(END_TIME/RECORD_STEP + 1, 1);
sum_REC_A2    = zeros(END_TIME/RECORD_STEP + 1, 1);
sum_REC_Pz    = zeros(END_TIME/RECORD_STEP + 1, 1);
sum_REC_PzA2  = zeros(END_TIME/RECORD_STEP + 1, 1);
sum_REC_PzA4  = zeros(END_TIME/RECORD_STEP + 1, 1);
sum_REC_mRNAz = zeros(END_TIME/RECORD_STEP + 1, 1);

sum_REC_l1    = zeros(END_TIME/RECORD_STEP + 1, 1);


% rate constants
k_xlate      = 0.1   /60;
k_xcribe     = 0.2   /60;
k_dec_mrna_a = 2.0   /60;
k_dec_mrna_z = 2.0   /60;
k_dec_a      = 0.05  /60;
k_dec_a2     = 0.05  /60;
k_dimer_a    = 1.0   /60;
k_sngl_a     = 0.01  /60;
k_rprs_a2    = 1.0   /60;
k_rprs_a4    = 0.01  /60;
k_dissoc_a2  = 10.0  /60;
k_dissoc_a4  = 0.2   /60;

k_l1         = 100.0 / 60;


% Setup Time Record Field
index = 1;
for i = 0 : RECORD_STEP : END_TIME,
    REC_time(index) = i;
    index = index + 1;
end    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastic Simulation Core Loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1 : REGRESSION
    
    % setup initial species values in the cell
%     mRNAa = 0;
    l1 = 0;
    
    A     = 0;
    A2    = 0;
    Pz    = PLASMID_COPIES;
    PzA2  = 0;
    PzA4  = 0;
    mRNAz = 0;
    rRNA  = 16; 
    RNAp  = 6; 

    % initialize internal variables
    sum_avals     = 0.0;
    index         = 1;
    curr_time     = 0.0;
  
    % main loop to simulate the system of reactions
    while( curr_time < END_TIME )
        
        % mRNAa Input Pattern (user customizable)
%         if( curr_time < 10800 )
%             mRNAa = 0;
%         elseif (curr_time < 21600 )
%             mRNAa = 30;
%         else
%             mRNAa = 0;
%         end    
        if( curr_time < 10800 )
            l1 = 0;
        elseif (curr_time < 21600 )
            l1 = 1;
        else
            l1 = 0;
        end    



        % random numbers used to determine which next reaction and when
        r1 = rand;    % when
        r2 = rand;    % which
                
        if( (r1 ~= 0) & (r2 ~= 0) )
            
            sum_avals = Calculate_Sum_Of_Avals;
            curr_time = curr_time + log(1/r1) / sum_avals;  % next time
            r2 = r2 * sum_avals;
            
            if( (REC_time(index) < curr_time) & (curr_time < END_TIME) )
                index = index + 1;
            end    
            tmp_sum = 0.0;
            eq_index = 0;
            
            % determine which reaction is next
            while( tmp_sum < r2 )
                tmp_sum = tmp_sum + AVALS( eq_index+1 );
                eq_index = eq_index + 1;    
            end    
            
            % execute reaction (i.e. update species values)
            EQ_Get_Next_States(eq_index);

            % record information about species values at this time point            
            REC_mRNAa(index) = mRNAa;
            REC_A(index)     = A;
            REC_A2(index)    = A2; 
            REC_Pz(index)    = Pz;
            REC_PzA2(index)  = PzA2;
            REC_PzA4(index)  = PzA4;
            REC_mRNAz(index) = mRNAz;
            REC_l1(index)    = l1;
           
        end
    end
    
    % Summarize the record field at each simulation
    index = 1;
    for i = 0 : RECORD_STEP : END_TIME,
        sum_REC_mRNAa(index) = sum_REC_mRNAa(index) + REC_mRNAa(index);
        sum_REC_A(index)     = sum_REC_A(index)     + REC_A(index);
        sum_REC_A2(index)    = sum_REC_A2(index)    + REC_A2(index); 
        sum_REC_Pz(index)    = sum_REC_Pz(index)    + REC_Pz(index);
        sum_REC_PzA2(index)  = sum_REC_PzA2(index)  + REC_PzA2(index);
        sum_REC_PzA4(index)  = sum_REC_PzA4(index)  + REC_PzA4(index);
        sum_REC_mRNAz(index) = sum_REC_mRNAz(index) + REC_mRNAz(index);
        sum_REC_l1(index) = sum_REC_l1(index) + REC_l1(index);
        index = index + 1;
    end    

end   
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting Figures            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
subplot(8,1,8);
plot(REC_time/60,sum_REC_l1/ REGRESSION);
hold on;
legend('l1');
hold off;

subplot(8,1,1);
plot(REC_time/60,sum_REC_mRNAa / REGRESSION);
hold on;
legend('mRNAa');
hold off;

subplot(8,1,2);
plot(REC_time/60,sum_REC_A / REGRESSION);
hold on;
legend('A');
hold off;

subplot(8,1,3);
plot(REC_time/60,sum_REC_A2 / REGRESSION);
hold on;
legend('A2');
hold off;

subplot(8,1,4);
plot(REC_time/60,sum_REC_Pz / REGRESSION);
hold on;
legend('Pz');
hold off;

subplot(8,1,5);
plot(REC_time/60,sum_REC_PzA2 / REGRESSION);
hold on;
legend('PzA2');
hold off;

subplot(8,1,6);
plot(REC_time/60,sum_REC_PzA4 / REGRESSION);
hold on;
legend('PzA4');
hold off;

subplot(8,1,7);
plot(REC_time/60,sum_REC_mRNAz / REGRESSION);
hold on;
legend('mRNAz');
hold off;