
function [tt,nn,aa] = gillespie(pars)

% Supplementary Matlab files for the paper:
%
% "Modeling and Simulation of IntraCellular Dynamics: Choosing an
% Appropriate Framework"
% by O.Wolkenhauer, M.Ullah, W.Kolch, K.-H.Cho
%
% Script file, comparing Gillespie's stochastic simulation (direct method)
% with the solution of the generalized mass action (GMA) model.
%
% We would appreciate a citation if these files are used by others.

n0 = pars{7};
nu = pars{9};
tf = pars{10};

rand('state',sum(100*clock)) % reset random number generator
t = 0;        % intialize time to zero
n = n0;        % initialize state     
tt = 0;
nn = n0;
aa = propensity(n0,pars);
while   t < tf
    a = propensity(n,pars); % compute propensities for all reaction channels
    astr = sum(a);          % sum propensities over reaction channels
    if ~astr
        t = tf;
    else
        tau = exprnd(inv(astr));       % generate 'tau' = time to next reaction 
        mu = min( find( cumsum(a) > astr*rand ) ); %generate 'mu' specifying next reaction Ru
        t = t + tau;                % update t 
        n = n + nu(mu,:);             % update n
    end       
    tt = [ tt; t ];
    nn = [ nn; n ];
    aa = [ aa; a ];
end
tt(end) = tf;  
nn(end) = nn( end - 1 );    
aa(end) = aa( end - 1 );

% computation of propensity(a)
function a = propensity(n,pars)
c = pars{3};
M = pars{4};
p = pars{5};
l = pars{6};
for mu=1:M
    if n(p{mu})
        h(mu) = prod(matnchoosek(n(p{mu}),l{mu}));  
    else
        h(mu) = 0;  
    end
end
a = c.*h;   % propensities