function cells = stoch_gma(pars)
%global n0 runs
% 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.

% Run the scriptfile which defines the model and simulation parameters.

%reactiondata; 

%rand('state',sum(100*clock)) % reset random number generator
n0 = pars{7};
runs = pars{end};

t_gilsp =  [];
for r=1:runs
    [t{r},n{r},a{r}] = gillespie(pars); % t: reaction times, n: populations
    astr{r} = sum(a{r},2);
    t_gilsp = union(t_gilsp,t{r});
    tau{r} = diff(t{r});
end

Lt = length(t_gilsp);
ensem_n = repmat(n0,[Lt,1,runs]);
ensem_a = repmat([0 0],[Lt,1,runs]);
for i=1:Lt
    for r=1:runs
        ind=find( t{r} <= t_gilsp(i) );
        ensem_n(i,:,r) = n{r}(ind(end),:); 
        ensem_a(i,:,r) = a{r}(ind(end),:); 
    end      
end

mean_n = mean(ensem_n,3); % mean of state over realisations
std_n = std(ensem_n,0,3); % standard deviation of state over realisations
mean_a = mean(ensem_a,3); % mean of state over realisations
std_a = std(ensem_a,0,3); % standard deviation of state over realisations

tspan = minmax(t_gilsp');
[t_gma,n_gma] = gma(tspan,pars);  % Solution of GMA

cells = {t,n,t_gma,n_gma,mean_n,std_n,t_gilsp,a,astr,mean_a};