% 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.

clear all,clc,close all

V = 1e-3;                  % volume of the compartment  (nL)
NAV = V*1e-9*6.02214199e23; %NAV = Avogadro's number times volume;
alpha = 1; beta = 1; gamma = 1;  % stoichiometric coefficients alpha beta gamma
p = {[1 2],3};        % participant(reactant) specie indicator
l = {[1 alpha],beta};        % stoichometric coefficients of reactants
k = [0.5 0.2];      % rate constants [(nM.s)^-1  s^-1  s^-1]
M = length(k);                  % Number of reaction channels

% Calculation of k' and c:

for u=1:M
    K(u) = sum(l{u});    % Ru reactant molecularity 
    kp(u) = k(u)/NAV^(K(u)-1); % 'particle' rate constant of Ru
    kp(u) = kp(u)*(1e9)^(K(u)-1);  % normalisation of kp(u) to have correct units
    c(u) = kp(u)*prod(matfactorial(l{u})); % stochastic rate constant
end 
% kp stands for kprime in the paper

% Step changes in all species in all reactions
%     S1     S2    S3     S4
nu = [ -1 -alpha  beta      0    % R1
       0  alpha -beta  gamma ]; % R2
n0 = [10 10 0 0];   % initial population 
S0 = 1e9*n0/NAV;       % initial concentration (nM)  
%S0 = [0.02 0.02 0 0];      % initial concentration (nM) 
%n0 = round(NAV*S0*1e-9);      % initial population 
tf = 400; % simulation end in seconds
runs = 50;  % number of realizations(simulation runs), over which to average