%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastically simulate enzyme functionality
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

%% Define the reaction

% S + E -> SE
% SE -> S + E
% SE -> P + E

%      S  E SE  P
lhs = [1, 1, 0, 0;
       0, 0, 1, 0;
       0, 0, 1, 0];
   
%      S  E SE  P
rhs = [0, 0, 1, 0;
       1, 1, 0, 0;
       0, 1, 0, 1];
   
c = [1.66e-3, 1e-4, 0.1]';
y0 = [301, 120, 0, 0]';

T_MAX = 50;
T_STEP = 0.01;

%% Time evolution plot

sim_log = ssas_gillespie(lhs, rhs, c, y0, T_MAX, T_STEP);

figure(1);

plot(sim_log(:,1), sim_log(:,2:end), '-');
legend('[S]','[E]', '[SE]', '[P]');
xlim([0 50]); ylim([0 300]);
xlabel('Time');
ylabel('# of molecules');

%% Run multiple trials

NUM_TRIALS = 100; % Number of trials. Reduce to decrease running time. Increase to obtain smoother plots

logs = zeros((T_MAX/T_STEP)+1, NUM_TRIALS);
for i = 1:NUM_TRIALS
    sim_log = ssas_gillespie(lhs, rhs, c, y0, T_MAX, 0.01);
    logs(:,i) = interp1q(sim_log(:,1), sim_log(:,2), (0:T_STEP:T_MAX)');
    
    if mod(i,10) == 0
        disp(i)
    end
        
    % Plot first 20 trials
    if i <= 20
        figure(2);
        plot(sim_log(:,1), sim_log(:,2), '-');
        hold on;
    end
end

legend('[P]');
xlabel('Time');
ylabel('# of molecules');

%% Mean/variance display

sim_avg = mean(logs,2);
sim_std = std(logs,0,2);
plot(0:T_STEP:T_MAX, sim_avg, '-', 'Color', 'red', 'LineWidth', 2);
plot(0:T_STEP:T_MAX, sim_avg+3*sim_std, '--', 'Color', 'green', 'LineWidth', 2);
plot(0:T_STEP:T_MAX, sim_avg-3*sim_std, '--', 'Color', 'green', 'LineWidth', 2);

%% PMF display

figure(3);

t = 10;
hist(logs((t/T_STEP)+1,:),0:301);
xlim([100 200])