# Lotka-Volterra simulation using Gillespie algorithm

22 vues (au cours des 30 derniers jours)
Serhat Unal le 22 Sep 2020
Commenté : Harley Day le 24 Sep 2020
&&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 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Harley Day le 24 Sep 2020
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]);
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

Connectez-vous pour commenter.

### Réponses (1)

Alan Stevens le 23 Sep 2020
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 commentairesAfficher AucuneMasquer Aucune
Serhat Unal le 23 Sep 2020
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:
and here they have used Y1 and Y2, but in my case I have several values for each for different timestep.
Harley Day le 24 Sep 2020
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

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Interactive Control and Callbacks dans Help Center et File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by