%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastically simulate autoregulatory gene network
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

%% Define the reaction

% G + P2 -> P2G
% P2G -> G + P2
% G -> G + R
% R -> R + P
% 2P -> P2
% P2 -> 2P
% R -> 0
% P -> 0

%      G P2 P2G R  P
lhs = [1, 1, 0, 0, 0;
       0, 0, 1, 0, 0;
       1, 0, 0, 0, 0;
       0, 0, 0, 1, 0;
       0, 0, 0, 0, 2;
       0, 1, 0, 0, 0;
       0, 0, 0, 1, 0;
       0, 0, 0, 0, 1];
   
%      G P2 P2G R  P
rhs = [0, 0, 1, 0, 0;
       1, 1, 0, 0, 0;
       1, 0, 0, 1, 0;
       0, 0, 0, 1, 1;
       0, 1, 0, 0, 0;
       0, 0, 0, 0, 2;
       0, 0, 0, 0, 0;
       0, 0, 0, 0, 0];
   
c = [1, 10, 0.01, 10, 1, 1, 0.1, 0.01]';
y0 = [10, 0, 0, 0, 0]';

T_MAX = 500; %Running time is 48 secs to simulate 100secs!, 468 for 500, 12802 for 10000
T_STEP = 0.01;

%% Sim
tic
sim_log = ssas_gillespie(lhs, rhs, c, y0, T_MAX, T_STEP);
toc
%% Time evolution plot
figure(1);

subplot(3,1,1);
plot(sim_log(:,1), sim_log(:,5), '-r');
ylabel('mRNA (R)');
ylim([0 2]);
subplot(3,1,2);
plot(sim_log(:,1), sim_log(:,6), '-g');
ylabel('Protein (P)');
subplot(3,1,3);
plot(sim_log(:,1), sim_log(:,3), '-b');
ylabel('Protein Dimer (P2)');

xlabel('Time');
samexaxis('xmt','on','ytac','join','yld',1)

%% PMF display

figure(2);

NUM_TRIALS = T_MAX/T_STEP;
[N,X]=hist(sim_log(:,3),0:60);
bar(X, N/NUM_TRIALS);
xlim([0 60]);
ylim([0 0.07]);
xlabel('# of Protein (P) Molecules');
ylabel('Probability');