% .. code2.m
% .. two-step gene expression model
% .. Noise in gene expression lectures

clear all
clc

% .. model parameters
kmRNA = 10;         % .. transcription rate; [] = #molecu/time
dmRNA = 10;         % .. mRNA degradation rate; [] = 1/time
kProt = 50;         % .. translation rate; [] = 1/time
dProt = 1;          % .. P degradation rate; [] = 1/time
tccyc = 30;         % .. cell cycle time; [] = minutes

% .. Stochastic dynamics simulation by Gillespie's algorithm
% .. initial values 
tmax = 10*tccyc;
tini = 0;
Ps = [0,0];         % .. Ps = [mRNA, Prot];
% .. optimizing code
tt = 100000;
Pstochastic = zeros(tt,2);
Pstochastic(1,:) = Ps;
ts = zeros(tt,1);
ts(1) = tini;
t = tini;

% ..
ss = 1;

% ..
while t < tmax
    
    % .. propensities
    a = [kmRNA, dmRNA*Ps(1), kProt*Ps(1), dProt*Ps(2)];
    a0 = sum(a);

    % .. determine time of next reaction
    r1 = rand;
    tau = -log(r1)/a0;
    t = t + tau;
    disp([t/tccyc, tmax/tccyc])
    
    % .. determine nature of next reaction
    r2 = rand;
    acumsum = cumsum(a)/a0;
    chosen_reaction = find(r2 <= acumsum,1);

    if chosen_reaction == 1;
        Ps(1) = Ps(1) + 1;
    elseif chosen_reaction == 2;
        Ps(1) = Ps(1) - 1;
    elseif chosen_reaction == 3;
        Ps(2) = Ps(2) + 1;
    else
        Ps(2) = Ps(2) - 1;
    end

    ss = ss + 1;
    ts(ss) = t;
    Pstochastic(ss,:) = Ps;
end

tss = (ts > 0);
ts = [tini;ts(tss)];
Pstochastic = [[0,0];Pstochastic(tss,:)];

% .. deterministic equations
Pd0 = [0,0];
options = [];
[t Pd] = ode23(@code2equations,ts,Pd0,options,kmRNA,dmRNA,kProt,dProt);

% .. plots
subplot(211)
plot(t/tccyc,Pstochastic,t/tccyc,Pd,'LineWidth',1.7)
axis([0,tmax/tccyc, 0, max(max(Pstochastic))])
xlabel('time/tcellcycle')
legend('det-mRNA','det-Prot','sto-mRNA','sto-Prot',0)
grid

% ..
subplot(212)
Pst = Pstochastic(:,2);
bar(min(Pst):max(Pst),histc(Pst,min(Pst):max(Pst)))
title(strcat('Deterministic value:', num2str(kmRNA*kProt/dmRNA/dProt)))

mean_st = mean(Pst)
CV = std(Pst)/mean(Pst)
Fano = var(Pst)/mean(Pst)
