Lotka-Volterra simulation using Gillespie algorithm

6 次查看（过去 30 天）
Serhat Unal
Serhat Unal 2020-9-22
评论： Harley Day 2020-9-24
&&calling this function in the script
function dx = myalg(t,x)
dx = zeros(3,1);
a = 0.25;
b = 0.0025;
d = 0.125;
dx(1) = 0; %food
dx(2) = a*x(1)*x(2)-b*x(2)*x(3); %R
dx(3) = b*x(2)*x(3)-d*x(3); %F
end
&&script begins here
clear all
clc
a = 0.25;
b = 0.0025;
d = 0.125;
tplot = zeros(1);
n = 0;
n_max = 1000000;
t = 0;
t_max = 30;
while t<t_max
[t x] = ode45(@myalg,[t t_max],[1 50 50]);
Rplot = x(:,2);
Fplot = x(:,3);
h = [a.*x(:,2),b.*x(:,3).*x(:,2),d.*x(:,3)];
ktot = sum(h);
r1 = randi([0, 4294967295]) / 4294967296.0;
tau = -(1./ktot)*log(r1);
r2 = randi([0, 4294967295]) / 4294967296.0;
mu = sum(r2*ktot <= cumsum(h));
t = t + tau;
switch mu
case 3
x(:,2) = x(:,2) + 1;
case 2
x(:,2) = x(:,2) - 1;
x(:,3) = x(:,3) + 1;
case 1
x(:,3) = x(:,3) - 1;
end
n = n + 1;
Rplot(n+1,:) = x(:,2);
Fplot(n+1,:) = x(:,3);
tplot(n+1) = t;
end
figure
hold on
stairs(tplot, Rplot)
stairs(tplot, Fplot)
legend({'$R$','$F$'}, 'interpreter','latex');
xlabel('time (minutes)', 'interpreter', 'latex');
ylabel('numbers of $R$ and $F$ molecules', 'interpreter', 'latex');
hold off
figure
plot(Rplot,Fplot);
xlabel('number of $R$ molecules', 'interpreter', 'latex');
ylabel('number of $F$ molecules', 'interpreter', 'latex');
%%I wonder why my code is not working, hope somebody can help me out.
I dont know how to recall the function, in other words where I should put [t x] = ode45(@myalg,[t t_max],[1 50 50])
and also why swith mu not works. Hope for help, thanks!
  1 个评论
Harley Day
Harley Day 2020-9-24
Hi Serhat,
I responded to your question on the other page. I hope my answer helps you. As for your code above, you seem to be trying to call a detministic simulation within the Gillespie algorithm's while loop. The Gillespie algorithm does not require any ODE solving to be carried out, so you should first remove the line
[t x] = ode45(@myalg,[t t_max],[1 50 50]);
from your code.
I'd love to be able to help more, but I can't work out what you're trying to do here. If you're just trying to simulate the system, please see my answer on the other page, as I've written a script that does just that.
I hope this helps. Feel free to respond if you need any more guidance.
Harley

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Alan Stevens
Alan Stevens 2020-9-23
You have some mismatched vector/matrix sizes. The following works, though you will need to check carefully to see if it does exactly what you want.
a = 0.25;
b = 0.0025;
d = 0.125;
n = 0;
n_max = 1000000;
t = 0;
t_max = 30;
[t, x] = ode45(@myalg,[t t_max],[1 50 50]);
h = [a.*x(:,2),b.*x(:,3).*x(:,2),d.*x(:,3)];
ktot = sum(h);
r1 = randi([0, 4294967295]) / 4294967296.0;
tau = -(1./ktot)*log(r1);
r2 = randi([0, 4294967295]) / 4294967296.0;
mu = sum(r2*ktot <= cumsum(h));
t = t + tau(:,1);
if mu==3
x(:,2) = x(:,2) + 1;
elseif mu==2
x(:,2) = x(:,2) - 1;
x(:,3) = x(:,3) + 1;
else
x(:,3) = x(:,3) - 1;
end
n = n + 1;
Rplot = x(:,2);
Fplot = x(:,3);
tplot = t;
figure
hold on
stairs(tplot, Rplot)
stairs(tplot, Fplot)
legend({'$R$','$F$'}, 'interpreter','latex');
xlabel('time (minutes)', 'interpreter', 'latex');
ylabel('numbers of $R$ and $F$ molecules', 'interpreter', 'latex');
hold off
figure
plot(Rplot,Fplot);
xlabel('number of $R$ molecules', 'interpreter', 'latex');
ylabel('number of $F$ molecules', 'interpreter', 'latex');
function dx = myalg(~,x)
dx = zeros(3,1);
a = 0.25;
b = 0.0025;
d = 0.125;
dx(1) = 0; %food
dx(2) = a*x(1)*x(2)-b*x(2)*x(3); %R
dx(3) = b*x(2)*x(3)-d*x(3); %F
end
  2 个评论
Serhat Unal
Serhat Unal 2020-9-23
I don't know if the code here does take into account iteration based on time. Because
we want to loop it for every timestep how every variable does change, but I dont think this code does that.
I can see that the while loop is removed and maybe should be there but calling the function should maybe be
outside the loop, I am not sure but I used this method from this link:
https://se.mathworks.com/matlabcentral/answers/224604-implementing-gillespie-s-algorithm
and here they have used Y1 and Y2, but in my case I have several values for each for different timestep.
Harley Day
Harley Day 2020-9-24
You've got to separate the ODE45 bits from the while loop of the Gillespie algorithm. They should not overlap at all.
The Gillespie algorithm is designed to simulate the stochastic system, and does not use any ODE solving methods at all. The following script runs the Gillespie algorithm to simulate the stochastic system:
%clear
%% Lotka reactions
%This gives figure 6 of Gillespie's 1977 paper.
M = 3; % Number of reaction pathways
N = 2; % Number of molecular species considered
%% STEP 0
% Initial number of molecules of species X(1) and X(2). In general, a
% vector of length N of initial population numbers for each species.
%--------------------------
% equation | rate
%--------------------|-----
% Xbar + Y1 -> 2*Y1 | c1
% Y1 + Y2 -> 2*Y2 | c2
% Y2 -> Z | c3
%--------------------------
Y1 = 1000;
Y2 = 1000;
X = 10;
% In general, a vector of length M of stochastic reaction constants.
c = [10/X,0.01,10]; % units: (expected number of) reactions per minute
Y1plot = Y1; % For plotting.
Y2plot = Y2;
t = 0; % Time variable.
t_max = 30;
tplot = zeros(1); % For plotting.
n = 0; % Reaction counter.
n_max = 1000000;
%% STEP 1
while t < t_max
% Number of molecular reactant combinations available in current state.
% In general, h is a vector of length M consisting of combinatorial
% functions of molecular population numbers X (which is of length N).
h = [X*Y1,Y1*Y2,Y2];
% a is the propensity of the reaction pathway in current state. In
% general, a vector of length M
a = h.*c;
% a0 is total propensity that anything happens. This number emerges
% more out of mathematical necessity than physical intuition.
a0 = sum(a);
%% STEP 2
r = rand(2,1);
tau = -log(r(1))/a0;
% Note 0<=mu<=M. This decides which reaction occurs. If mu=0, no
% reaction occurs. The switch statement does nothing in this case.
mu = sum(r(2)*a0 <= cumsum(a));
%% STEP 3
t = t + tau;
% Adjust population levels based on reaction formula(s).
switch mu
case 3
Y1 = Y1 + 1;
case 2
Y1 = Y1 - 1;
Y2 = Y2 + 1;
case 1
Y2 = Y2 - 1;
end
n = n + 1;
% At this point, all the physics has been simulated, it only remains to
% record the new values of t and X to vectors to we can plot it later.
Y1plot(n+1,:) = Y1;
Y2plot(n+1,:) = Y2;
tplot(n+1) = t;
end
figure
hold on
stairs(tplot, Y1plot)
stairs(tplot, Y2plot)
legend({'$Y_1$','$Y_2$'}, 'interpreter','latex');
xlabel('time (minutes)', 'interpreter', 'latex');
ylabel('numbers of $Y_1$ and $Y_2$ molecules', 'interpreter', 'latex');
hold off
figure
plot(Y1plot,Y2plot);
xlabel('number of $Y_1$ molecules', 'interpreter', 'latex');
ylabel('number of $Y_2$ molecules', 'interpreter', 'latex');
Note that we don't call ode45 at all in the above script.
If you want to simulate the deterministic system, you just give the ODEs to the ode45 function to solve numerically. This is done in the following script:
[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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by