function varargout = stoch(varargin)
%
% STOCH stochastic simulation of a biochemical network
% [Ts,Ns] = STOCH(X0,N0,K,C,D,L,TF) performs stochastic simulation a
% biochemical network composed of M elementary reactions, using Gillespie
% algorithm. The input and output arguments have the following meaning:
%
% X0: Column vector of initial concentrations of all the species involved
% N0: Column Vector of initial populations of all the species involved
% K: Row vector of rate constants of all elementary reactions
% C: Row vector of stochastic rate constants of all elementary reactions
% D: Stoichiometry matrix with rows correspoding to species and columns
% corresponding to reaction channels
% L: Stoichiometry matrix for reactants only such that L = -D.*(D < 0);
% TF: Final time of simulation
% Ts: Row vector of time points of reaction events
% Ns: Matrix of output concentrations with a column for each time point.
%
% [Ts,Ns] = STOCH(X0,N0,K,C,D,L,TF,R) performs R runs of stochastic
% simulation. The outputs Ts and Ns are cell arrays where each cell
% corresponds to one run.
%
% [Ts,Ns,TT,NBAR] = STOCH(X0,N0,K,C,D,L,TF,R) where R>1, also returns
% the ensemble of times in the row vector TT and the mean of Ns over R
% runs in matrix NBAR of the size of Ns.
%
% STOCH(X0,N0,K,C,D,L,TF,1) is the same as STOCH(X0,N0,K,C,D,L,TF).
% preliminary operations:
args = varargin;
if nargin<6, args{6} = 1; end % single run by default
[n0,c,d,l,tf,runs] = args{:}; % parse inputs
oneszk = ones(size(c));
i1 = l==1; % reactions of type: A->B, A+B->AB
i2 = l==2; % reactions of type: A+A->AA
stop = tf - eps(tf); % simulation stop time
nOut = nargout;
% stochastic part:
if runs < 2 % single run
[yout{1:2}] = gillespie; % run gillispie
else % multiple runs
[tg,ng] = deal(cell(1,runs));
for i = 1:runs, [tg{i},ng{i}] = gillespie; end % run gillispie
yout = {tg,ng};
if nOut>2
tt = unique([tg{:}]); % record times
numeltt = numel(tt); % record populations
nn = zeros(numel(n0), numeltt, runs);
for i = 1:numeltt
for j = 1:runs
id = find(tg{j} <= tt(i),1,'last');
nn(:,i,j) = ng{j}(:,id);
end
end
yout(3:4) = {tt,mean(nn,3)}; % append population mean
end
end
% Output:
varargout = yout;
% gillespie algorithm (direct method):
function [tt,nn] = gillespie
t = 0; % initial time
n = n0; % initial population
tt = [];
nn = [];
rand('state', sum(100*clock)); % reset random number generator
while t <= stop
tt = [tt t]; % record time
nn = [nn n]; % record population
m = n(:,oneszk); % replicate n : size(m,2) = size(k)
b = double(~l);
b(i1) = m(i1); % reactions of type: A->B, A+B->AB
b(i2) = m(i2).*(m(i2)-1)/2; % reactions of type: A+A->AA
a = c.*prod(b); % propensity
astr = sum(a);
if ~astr, break, end % substrate utilised
tau = -1/astr*log(rand); % time to next reaction
u = find(cumsum(a)>=astr*rand,1);% index of next reaction
n = n + d(:,u); % update population
t = t + tau; % update time
end
tt = [tt tf];
nn = [nn nn(:,end)];
end
end