How do I average multiple sets of data created by Gillespie's Algorithm?

Hello everyone,

I am pretty to new to compute science and only know some basic stuff such as the if statement. I tried to implement Gillespie's Algorithm and average all the data to get a plot. My function is as follows:

function [Gilm, Gilsm, Gilp, Gilsp] = GilUnreg

% simulate the unregulated gene expression
% -> mRNA (with production rate kr)
% -> protein (with production rate kp)
% mRNA -> (with degradation rate delr)
% protein -> (with degradation rate delp)

n = 800;

% set up data storage
t = zeros(1, n);
m = zeros(1, n); % number of mRNA
p = zeros(1, n); 
R = rand(1, n); % R determines which reaction should occur

% parameters
kr = 0.14; % kr = k0
kp = 0.14; % kp = k1
delr = log (2)/50; % delr = delm
delp = log (2)/50;

% initial conditions
m(1) = 0;
p(1) = 0;
t(1) = 0;

for i = 1:(n-1)
    
    % define propensities
    a1 = kr;
    a2 = kp*m(i);
    a3 = delr*m(i);
    a4 = delp*p(i);
    a = a1 + a2 + a3 + a4;
    
    % determine which reaction would happen
    if (R(i) >= 0) && (R(i) <= a1/a)
        m(i+1) = m(i) + 1;
        p(i+1) = p(i);
    elseif (R(i) > a1/a) && (R(i) <= (a1 +a2)/a)
        m(i+1) = m(i);
        p(i+1) = p(i) + 1;
    elseif (R(i) > (a1 + a2)/a) && (R(i) <= (a1 + a2 + a3)/a)
        m(i+1) = m(i) - 1;
        p(i+1) = p(i);
    elseif (R(i) > (a1 + a2 +a3)/a) && (R(i) <= 1)
        m(i+1) = m(i);
        p(i+1) = p(i) - 1;
    end
    
    % determine how long the reaction takes
    dt = -1/a*log(1-rand);
    t(i+1) = t(i) + dt;
    
end

end

If I run this code once, I will get a set of data with time, mRNA and protein. If I run this code again, I will get another set of data with time, mRNA and protein. I want to run it many times, and my question is then: how do I average all the data and get a graph of mRNA with respect to time and a graph of protein with respect to time?

[Mario Malic]

I don't see your output arguments defined in the code. This would create an array with 50 rows in which each one contains data from a function call. This would work if the output arguments are row vectors, which I can't check.

for ii = 1:1:50
    [Gilm(ii,:), Gilsm(ii,:), Gilp(ii,:), Gilsp(ii,:)] = GilUnreg
end

You can use function mean to get average of Gilm across all function calls and the others in a same way.

GilmAvg = mean (Gilm);

You can use plot function for plotting, you have to check that you have equal number of elements for your x and y axes. Values in your time vector vary with each run, so I am not completely sure how would you like to average that. Should you interpolate Gilm values of each function call for averaged time vector or do it some other way, it's up to you.

Note:

if (R(i) >= 0) && (R(i) <= a1/a) 
elseif (R(i) > a1/a) && (R(i) <= (a1 +a2)/a)
elseif (R(i) > (a1 + a2)/a) && (R(i) <= (a1 + a2 + a3)/a)
elseif (R(i) > (a1 + a2 +a3)/a) && (R(i) <= 1)
else % It would be good to include else block,
     % if all of the above conditions are not fullfilled
     % then it's not determined which reaction happens (if it happens)
     fprintf("Reaction is not determined for %f", R(i));
end

[Piyush Aggarwal]

For the implementaion of algorithm - The following file exchange link contains a code file with implementation of the algorithm.

https://www.mathworks.com/matlabcentral/fileexchange/34707-gillespie-stochastic-simulation-algorithm

This post also might turn to be helpful:

https://www.mathworks.com/matlabcentral/answers/224604-implementing-gillespie-s-algorithm