Asked by Argento
on 19 Jun 2015

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)

## 3 Comments

## Harley Day (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/224604-implementing-gillespie-s-algorithm#comment_296997

## Harley Day (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/224604-implementing-gillespie-s-algorithm#comment_670153