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

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

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

도움말 센터File Exchange에서 Chemistry에 대해 자세히 알아보기

제품

질문:

2020년 12월 24일

댓글:

2021년 1월 3일

Community Treasure Hunt

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

Start Hunting!

Translated by