% 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