MATLAB Answers

Implementing Gillespie's Algorithm.

293 views (last 30 days)
Argento
Argento on 19 Jun 2015
Edited: Harley Day on 1 Oct 2020
Hello everyone,
Being pretty new to Matlab, I've been struggling trying to implement Gillespie's Algorithm (1977). Truth be told, I am still somewhat confused by certain aspects of the algorithm itself (such as the calculation of the propensity function). I've been tying to stick pretty close to the methods outlined in his paper.
Below I posted my attempt at coding fig.6 from his (1977) paper, but I've been getting nowhere.
Could someone with some experience with the algorithm tell me what I'm doing wrong? Many thanks.
clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Species and reaction channels %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = 2;
N = 1;
%%%%%%%%%%%%%%
%%%STEP 0 %%%
%%%%%%%%%%%%%%
c1X = 5; %% ???? %%
c1 = 5; %% ???? %%
c2 = 0.005;
Y = 10;
X1 = Y;
t = 0; %%% Time %%%
n = 0; %%% Counter %%%
t_max = 100; %%% Max. number of iterations %%%
T = zeros(1, t_max + 1); %%% Time values array (for plotting purposes) %%%
XX = zeros(1, t_max + 1); %%% X values array %%%
Yc = zeros(1, t_max + 1); %%% Y values array (for plotting purposes) %%%
%%%%%%%%%%%%%%
%%%STEP 1 %%%
%%%%%%%%%%%%%%
while n < t_max
h1 = X1*Y;
h2 = Y*(Y-1)*0.5;
a = zeros(1, 2);
a1 = h1*c1;
a2 = h2*c2;
a(1, 1) = a1;
a(1, 2) = a2;
a0 = sum(a);
%%%%%%%%%%%%%%
%%%STEP 2 %%%
%%%%%%%%%%%%%%
r1 = rand();
r2 = rand();
tau = -log(r1)*(1/a0);
t = t + tau;
if r2*a0 <= sum(a);
Y = Y + h1;
Yc(1, n+1) = Y;
else
Y = Y - h2;
Yc(1, n+1) = Y;
end
%%%%%%%%%%%%%%
%%%STEP 3 %%%
%%%%%%%%%%%%%%
t = t + tau;
T(1, n+1) = t;
n = n + 1;
end
plot(T, Yc)
  4 Comments