% .. code1.m
% .. simplest gene expression model
% .. Noise in gene expression lectures

clear all
clc

% .. model parameters
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 = 30*tccyc;
tini = 0;
Ps = 0;         % .. Ps = [Prot];
% .. optimizing code
tt = 100000;
Pstochastic = zeros(tt,1);
Pstochastic(1,:) = Ps;
ts = zeros(tt,1);
ts(1) = tini;
t = tini;

% ..
ss = 1;

% ..
while t < tmax
    
    % .. propensities
    a = [kProt, dProt*Ps(1)];
    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;
    else
        Ps(1) = Ps(1) - 1;
    end

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

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

% .. deterministic equations
Pd0 = 0;
options = [];
[t Pd] = ode23(@code1equations,ts,Pd0,options,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-Prot','sto-Prot',0)
grid
% ..
subplot(212)
Pst = Pstochastic;
bar(min(Pst):max(Pst),histc(Pst,min(Pst):max(Pst)))
grid
title(strcat('Deterministic value:', num2str(kProt/dProt)))
xlabel('Proteins')
ylabel('Pn')

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