how to generate (multiple) realisations using Gillespie Algorithms

1 visualizzazione (ultimi 30 giorni)
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 Commenti
Star Strider
Star Strider il 3 Gen 2021
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.)
yamen alharbi
yamen alharbi il 3 Gen 2021
Thank you so much
Here is the first paper by Gillespie
and I have uploaded a copy of the same paper
Many Thanks again

Accedi per commentare.

Risposte (0)

Categorie

Scopri di più su Chemistry in Help Center e File Exchange

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by