%
% gillrunner.m   Steve Cox, October 6, 2006
%
% A scaled  down implementation of the Gillespie Simulation of 
% of McAdams and Arkin PNAS 1997 94 814-819
%
% Compute and display dimer stats for the multiple runs of the
% self-regulated gene.
%
% 8 Reactions:
%
% R+PR -> RPo   % RNAP promoter binding and opening 
% R+PR <- RPo   % RNAP promoter closing and unbinding
% RPo -> 10P    % mRNA production from RPo
% P + P -> D    % dimerization
% P + P <- D    % dedimerization
% P ->          % protein degradation
% PR + D -> I   % D steals PR
% I -> PR + D   % D falls off of PR
%
%  usage    gillrunner(tfin, tinc, nr, x, k)
%
%  where  	tfin = time at which to terminate
%
%               tinc = increment of time in revealed stats
%
%               nr = number of runs
%
%               x = initial `counts' of [R PR RPo P D I]
%
%               k = 8 propensities
%
%  e.g.,  gillrunner(10,.1,3,[10 1 0 0 0 0], [2 1 4 2 .5 0.05 .01 0.005])
%

function gillrunner(tfin,tinc,nr,x0,k)

tvec = 0:tinc:tfin;

% rtab{k}(1:2:end) = indicies of reaction species
% rtab{k}(2:2:end) = actions for species above

rtab = {[1 -1 2 -1 3 1]      % R+PR -> RPo
        [1 1 2 1 3 -1]       % RPo -> R + PR
        [1 1 2 1 3 -1 4 10]  % RPo -> 10P +R +PR
        [4 -2 5 1]           % P + P -> D
        [4  2 5 -1]          % D -> P + P
        [4 -1]               % P ->
        [2 -1 5 -1 6 1]      % PR + D -> I
        [2  1 5  1 6 -1]};   % I -> PR + D


for j=1:nr	% make nr Gillespie runs

    disp(['Proceeding with run ' num2str(j)])  % speak to user

    [t,x] = mygill(tfin, rtab, x0, k);        % record time and dimer count 

    X(j,:) = interp1(t,x,tvec);                % interpolate and store

end

m = mean(X,1);		% compute the ensemble (over runs not time!) mean

s = std(X,1);           % compute the ensemble mean standard deviation

plot(tvec,m)
hold on
plot(tvec,m+s,'r')
plot(tvec,m-s,'r')
hold off
xlabel('time','fontsize',16)
ylabel('dimer count','fontsize',16)

return

%
%  straight Gillespie on the self-regulating gene
%
%  usage  [t,dimer] = mygill(tfin, rtab, x, c)
%
%  where   tfin, rtab, x, and k are as described above
%          and t and dimer are the resulting time and dimer count
%

function [t,dimer] = mygill(tfin, rtab, x, c)

t = zeros(1,1e5);   % preallocate space for t and dimer, for speed
dimer = t;
t(1) = 0;
dimer(1) = x(5);

i = 1;

while t(i) < tfin,		% proceed until tfin

    h = update(x);              % # of distint reactant combinations
    a = c.*h;			% reaction probabilities
    a0 = sum(a);
    
    r1 = rand;
    i = i + 1;
    tau = log(1/r1)/a0;		% the time to next reaction
    t(i) = t(i-1) + tau;	% record time

    r2 = rand;
    b = cumsum(a)/a0;
    reac = find(b>r2,1);		% the reaction number 

    metabs = rtab{reac}(1:2:end);       % the associated reactants
    action = rtab{reac}(2:2:end);       % the associated action
    x(metabs) = x(metabs) + action;

    dimer(i) = x(5);		% record dimer count

end

disp(['    took ' num2str(i) ' Gillespie steps'])  % talk to user

t = t(1:i);		% only return the computed t and dimer vals
dimer = dimer(1:i);

return

%
% compute the number of distinct ways that each reaction may occur
%

function h = update(x)
h = zeros(1,8);
h(1) = x(1)*x(2);
h(2) = x(3);
h(3) = x(3);
h(4) = x(4)*(x(4)-1)/2;
h(5) = x(5);
h(6) = x(4);
h(7) = x(2)*x(5);
h(8) = x(6);