how to generate (multiple) realisations using Gillespie Algorithms
2 views (last 30 days)
Show older comments
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
Star Strider
on 3 Jan 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.)
Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!