how to generate (multiple) realisations using Gillespie Algorithms

Hello, I would like to simulate this reaction
a->X
X-> X
at rates c. here is a code was written by Desmond Higham and I tried to adopted it to this reaction.
I would like to simulate different realisations . Could you please shed a light on this.
Your help is highly aprrecaited.
clf
clear all;
clc
rand('state',100)
%stoichiometric matrix
V = [1 -1];
%%%%%%%%%% Parameters and Initial Conditions %%%%%%%%%
nA = 6.023e23; % Avagadro's number
vol = 1e-15; % volume of system
X = zeros(1,1);
X(1) = 100; % molecules of substrate
c(1) = .4; c(2) =.1; c(3)=8; c(4)=.3;
t = 0;
tfinal = 100;
count = 1;
tvals(1) = 0;
Xvals(:,1) = X;
while t < tfinal
a(1) = c(1);% calculate propensity functions
a(2) = c(2)*X(1);
asum = sum(a);
j = min(find(rand<cumsum(a/asum))); % select which reaction
tau = log(1/rand)/asum; % select time
X = X + V(:,j); %update X
count = count + 1;
t = t + tau; %Update time
tvals(count) = t;
Xvals(:,count) = X;
end
%%%%%%%%%%% Plots %%%%%%%%
L = length(tvals);
tnew = zeros(1,2*(L-1));
tnew(1:2:end-1) = tvals(2:end);
tnew(2:2:end) = tvals(2:end);
tnew = [tvals(1),tnew];
%
Svals = Xvals(1,:);
ynew = zeros(1,2*L-1);
ynew(1:2:end) = Svals;
ynew(2:2:end-1) = Svals(1:end-1);
plot(tnew,ynew,'g-')
hold on

5 Comments

Relevant information on the Gillespie algorithm would be helpful. I didn’t find anything that related to it that had sufficient detail when I did an Interweb search several days ago. (Obviously, I have no prior experience with it.)
Thanks for your replay
Here are some links
The Gillespie started here and it has the algoriyhm in Figure 2.
and here are some useful
Thanks again
The only one that appears to be accessable is the SIAM paper. I’ll download it now and read it later. (The first is behind a paywall, the second offers only short descriptions and no links to other information.)
Thank you so much
Here is the first paper by Gillespie
and I have uploaded a copy of the same paper
Many Thanks again

Sign in to comment.

Answers (0)

Categories

Find more on Chemistry in Help Center and File Exchange

Products

Asked:

on 24 Dec 2020

Commented:

on 3 Jan 2021

Community Treasure Hunt

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

Start Hunting!