Lotka-Volterra simulation using Gillespie algorithm

&&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!

[Harley Day]

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

Answer 1

[Alan Stevens]

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

[Serhat Unal]

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]

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