%
% smcregio.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 and contrast with the deterministic solution
%
% 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    smcregio(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.,  smcregio(10,.1,3,[10 1 0 0 0 0], [2 1 4 2 .5 0.05 .01 0.005])
%

function smcregio(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,xout] = gomcreg(tfin, rtab, x0, k);        % record time and dimer count 
    if nr == 1,
       tvec = t;
       X([1:6]+(j-1)*6,:) = xout;
    else
       tmp = interp1(t,xout',tvec);      % interpolate and store
       X([1:6]+(j-1)*6,:) = tmp';      % interpolate and store
    end

end

m = zeros(6,length(tvec));
s = m;

for j=1:6,
    m(j,:) = mean(X(j:6:end,:),1);  % compute ensemble mean
    s(j,:) = std(X(j:6:end,:),1);  % compute ensemble mean
end

[t,x] = ode23(@(t,x) geneode(t,x,k),[0 tfin],x0');

if nr == 1

   slab = {'R' 'PR' 'RPo' 'P' 'D' 'I'};

   for j=1:6
       subplot(2,3,j)
       plot(tvec,m(j,:),'.')
       hold on
       plot(t,x(:,j),'k')
       xlabel('t')
       ylabel(slab{j})
       ax = axis;
       ax(2) = tfin;
       axis(ax)
       hold off
   end

else

   plot(tvec,m(5,:))
   hold on
   plot(tvec,m(5,:)+s(5,:),'r')
   plot(tvec,m(5,:)-s(5,:),'r')
   plot(t,x(:,5),'k')
   xlabel('t')
   ylabel('dimer')
   legend('mean','mean+std','mean-std','ode','location','se')
   hold off

end

return

%  x = [R PR RPo P D I]

function dx = geneode(t,x,k)
dx(1,1) = -k(1)*x(1)*x(2) + k(2)*x(3) + k(3)*x(3);
dx(2,1) = -k(1)*x(1)*x(2) + k(2)*x(3) + k(3)*x(3) - k(7)*x(2)*x(5) + k(8)*x(6);
dx(3,1) = k(1)*x(1)*x(2) - k(2)*x(3) - k(3)*x(3);
dx(4,1) = 10*k(3)*x(3) - k(4)*x(4)^2 + 2*k(5)*x(5) - k(6)*x(4);
dx(5,1) = k(4)*x(4)^2/2 - k(5)*x(5) - k(7)*x(2)*x(5)  + k(8)*x(6);
dx(6,1) = k(7)*x(2)*x(5) - k(8)*x(6);

%
%  straight Gillespie on the self-regulating gene
%
%  usage  [t,dimer] = gomcreg(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,xout] = gomcreg(tfin, rtab, x, c)

t = zeros(1,1e5);   % preallocate space for t and dimer, for speed
xout = zeros(6,1e5); 
t(1) = 0;
xout(:,1) = x';

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;

    xout(:,i) = x';		% 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
xout = xout(:,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);