Implementing Gillespie's Algorithm.

41 views (last 30 days)
Argento
Argento on 19 Jun 2015
Answered: Jeremy Huard on 1 Jun 2022
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
David Goodmanson
David Goodmanson on 17 Feb 2019
Hi Harley,
quite a tour de force.
I have often felt that this site should allow voting for comments as well as answers, for the simple reason that comments can sometimes be better and more insightful than answers.
Serhat Unal
Serhat Unal on 23 Sep 2020
why is it Y1,Y2=1000,1000 and how would the code for Lotka reactions look like
if we solve the ODE's using ODE45 in a seperate m-file, the we would have
three columns and many rows for Y1,Y2 and X for each time-step.

Sign in to comment.

Answers (2)

Harley Day
Harley Day on 24 Sep 2020
Edited: Harley Day on 1 Oct 2020
Hi Serhat,
I assume you're talking about the Lotka reactions (equations 36 of the Gillespie paper).
We set
Y1 = 1000;
Y2 = 1000;
as our initial conditions. This choice of initial condition is completely arbitrary. You can change this is you like. In fact, if you copy paste the Lotka reaction code from above, and change the intial condition, you can see how/if it changes the dynamics.
The second part of your question seems to be about finding the deterministic equivalent of the Lotka example. The way we produce the detministic equivalent system is to look at each of the reactions in the system, converting each of these into terms in the ODE. I'll take us through the process for this example:
Let's just look at the first raction: . To have this reaction occur, we need the molecule to meet the molecule . In the deterministic system, we assume that the numbers of these molecules is represented by a continuous, single valued function (which remains constant so I'm going to drop the notation for ) and , which is a reasonable assumption when these numbers remain fairly large and the system is well-mixed.
The number of ways in which the molecules of can combine with the molecules of is , and we produce two molecules as a result (though one is just a reactant being released again so we only gain one molecule overall). This reaction occurs at rate . From all this information, we can write the first term of our ODE system:
We just have to work through the remaining reactions to fill in the remaining terms of the ODE. Note that the second reaction on the list will produce two ODE terms, since it uses up a molecule and produces a molecule.
When we've finished this process, we obtain the coupled system of ODEs:
These ODEs can be straighforwardly solved using ODE45 as you suggest. Here's a quick script to do just that, applying the initial conditions and , we discussed above.
[t, y] = ode45 ( @Lotka_deterministic, [0, 5], [1000, 1000] );
figure;
hold on;
plot ( t, y );
legend( {'$Y_1=1000$','$Y_2=1000$'}, 'interpreter', 'latex' );
xlabel('time (minutes)', 'interpreter', 'latex');
ylabel('mean number of $Y_n$ molecules', 'interpreter', 'latex');
xlim ( [0, 5] );
function dy = Lotka_deterministic(t, y)
Xbar = 10;
c_1 = 10/Xbar;
c_2 = 0.01;
c_3 = 10;
dy = [c_1*Xbar*y(1) - c_2*y(1)*y(2);
c_2*y(1)*y(2) - c_3*y(2)];
end
Now, clearly we're not seeing any dynamics here (note that the solution in blue is hidden behind the solution in red so we can't see it). This is because the system is at an (unstable) steady state, so the ODE45 algorithm evaluates the ODE at and discovers that and . You can confirm this using the code by running
Lotka_deterministic(0, [1000; 1000])
dy =
0
0
However, if you purturb the initial conditions by even a small amount, you will see limit cycle behaviours emerging. For example setting instead:
[t, y] = ode45 ( @Lotka_deterministic, [0, 5], [1000.1, 1000] ); % note the small change to the initial condition for Y_1
figure;
hold on;
plot ( t, y );
legend( {'$Y_1=1000$','$Y_2=1000$'}, 'interpreter', 'latex' );
xlabel('time (minutes)', 'interpreter', 'latex');
ylabel('mean number of $Y_n$ molecules', 'interpreter', 'latex');
xlim ( [0, 5] );
function dy = Lotka_deterministic(t, y)
Xbar = 10;
c_1 = 10/Xbar;
c_2 = 0.01;
c_3 = 10;
dy = [c_1*Xbar*y(1) - c_2*y(1)*y(2);
c_2*y(1)*y(2) - c_3*y(2)];
end
So why does this happen? Why does the stochastic simulation using the Gillespie algorithm show us limit-cycling behaviour even if we put the initial conditions on a steady state, while the deterministic equivalent just stays at steady state? The answer is that the random fluctuations in the numbers of molecules due to the random occurance of the reactions rapidly moves the system away from the steady state, pushing it towards the limit cycle.
The Gillespie paper uses this example specifically because it shows this counterintiutive behaviour: The deterministic system can show us a steady state, whereas the stochastic simulation shows us that the steady state is quickly destroyed by stochastic noise. The point Gillespie was trying to make was that the introduction of Stochastic noise can fundamentally change the behaviours of whatever system you're studying.
I hope this helps. Let me know if there are any points which are unclear. :)

Jeremy Huard
Jeremy Huard on 1 Jun 2022
While this post is about understanding Gillespie's algorithm, it might be interesting for anyone landing on this page to know that this algorithm is one of the stochastic solvers implemented in SimBiology.
Figure 6 of Gillespie's paper can be reproduced with the following code.
Please note that SimBiology requires deterministic rate constants and will scale them automatically to compute stochastic rate constants.
% create SimBiology model
model = sbiomodel('reaction29');
c = model.addcompartment('c',1);
% species
c.addspecies('Xbar',5,Constant=true);
y = c.addspecies('Y',3000);
c.addspecies('Z',0);
% stochastic rate constants
c1 = 1;
c2 = 0.005;
% deterministic rate constants
volume = c.Value;
model.addparameter('k1',volume*c1);
model.addparameter('k2',volume*c2/2);
% reactions
r1 = model.addreaction('Xbar + Y -> 2 Y');
kin1 = r1.addkineticlaw('MassAction');
kin1.setparameter( 'Forward Rate Parameter','k1');
r2 = model.addreaction('2 Y -> Z');
kin2 = r2.addkineticlaw('MassAction');
kin2.setparameter( 'Forward Rate Parameter','k2');
verify(model);
% run stochastic simulations
cs = model.getconfigset();
cs.StopTime = 5;
cs.RuntimeOptions.StatesToLog = 'Y';
cs.SolverType = 'ssa';
cs.SolverOptions.LogDecimation = 10; % 10 rpd plot
y.Value = 3000;
results1 = sbiosimulate(model, cs);
y.Value = 10;
results2 = sbiosimulate(model, cs);
% run deterministic simulations
cs.SolverType = 'ode15s';
y.Value = 3000;
results1d = sbiosimulate(model, cs);
y.Value = 10;
results2d = sbiosimulate(model, cs);
% plot results
colors = colororder;
figure;
hold on;
stairs(results1.Time, results1.Data, Color=colors(1,:), LineStyle='--');
plot(results1d.Time, results1d.Data, Color=colors(1,:), LineStyle='-');
stairs(results2.Time, results2.Data, Color=colors(2,:), LineStyle='--');
plot(results2d.Time, results2d.Data, Color=colors(2,:), LineStyle='-');
yline(1000);
legend("Y0 = 3000, stochastic","Y0 = 3000, deterministic",...
"Y0 = 10, stochastic", "Y0 = 10, deterministic")
ylabel("number of Y molecules")
xlabel("time")
box off;
set(gca,"XLimitMethod","padded","YLimitMethod","padded");

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!