Currently I am doing some Teaching Assistance for an advanced Math Biology / Biophysics class. Wanting to point our students towards a simple example of a Gillespie algorithm (check references by Gillespie and textbooks for a description of the algorithm)  in MatLab, I found that google served up one fully developed set of boxed functions (not good for understanding how this really works), or page-long tutorials (homework is due in a week or so – not good either). As we all like copy&paste style code snippets, here is … a code snippet that simulates the simplest reaction scheme I could come up with on the spot. $A \leftrightharpoons B$, with forward reaction ($A\rightarrow B$) rate $k_+$ and reverse reaction ($A\leftarrow B$) rate $k_-$. Additionally a constant number of particles, $A+B=N$, was assumed, allowing the replacement $B = N-A$, which reduces the system to one dimension.

Was this code helpful for you? How do you implement your Gillespie style simulations? Did you know that the Gillespie algorithm was “invented” by several people independently, and is also called Kinetic Monte Carlo, Doob-Gillespie, or Lebowitz algorithm?

(Also sorry that wordpress ate my indentation. Just copy the code into your MatLab editor, hit CTRL+A, then CTRL+I on a PC, or COMMAND+A, COMMAND+I on a Mac and it’s fixed.)

% This simulates low particle numbers going through the following reaction
% scheme:
%
% A -> B with rate k_plus, propensity A*k_plus
% B -> A with rate k_minus, propensity B*k_minus = (N-A)*k_minus
%
% The particle number is NN, and summing the particles in A and B up gives
% the constant number NN. This means no particles get destroyed or created,
% their number is conserved. Only conversion between the two states occurs.
%
% We also assume that the reaction rate to go A->B (k_plu) and from B<-A
% (k_minus) are constants.

%%% Model parameters and initial conditions

% Model parameters, i. e. the constant reaction rates
k_plus = 10;
k_minus = 1;

% state changes for the two reactions
delta_AA = [-1,+1];

% Initial conditions
AA = 100; BB = 100; tt = 0;
% Determien the total number of partiles
NN = AA + BB;
% We really only have to follow AA, because BB = NN - AA

%%% Set the parameters for the execution of the simulation
max_steps = 10000; % maximal number of reactions to simulate
max_time = 1; % maximal time the simulation should run for

% Preallocate the array that will contain the AA counts
AA_array = zeros(1,max_steps+1);
AA_array(1) = AA; % Apply initial conditions to first element
% Preallocate the array that will contain the time points
tt_array = zeros(1,max_steps+1);

% Predraw the random numbers
rand_nums = rand(2,max_steps);
% One row for the reaction times, one for the reaction channel

%%% Actual simulation
%initialize the reaction counter
r_counter = 0;
while ( r_counter < max_steps && tt < max_time )

r_counter = r_counter + 1; % Increase reaction counter by 1

% Calculate the current propensities
aa = [AA_array(r_counter).*k_plus, ...
(NN-AA_array(r_counter)).*k_minus];

% get indices of nonzero propensities
valid_inds = aa>0;
% Get valid aa and reaction channels
valid_aa = aa(valid_inds);
valid_delta_AA = delta_AA(valid_inds);
% Sum of all propensities
aa_sum = sum(valid_aa);

% Calculate random waiting taime
delta_t = -log(rand_nums(1,r_counter))./aa_sum;

% Construct interval to pick next reaction
aa_interval = cumsum(valid_aa); % Cumulative sum
aa_interval = aa_interval./aa_sum; % normalize into range [0,1]

% Find a random index using the weights contained in the normalized
% interval. This gives the first index out of the valid indices that
% exceeds the random number for this step.
next_reaction_ind = find(aa_interval>rand_nums(2,r_counter), ...
1,'first');

% Execute the reaction that is associated with this index
chosen_delta_AA = valid_delta_AA(next_reaction_ind);

% Apply the time and AA changes
tt_array(r_counter+1) = tt_array(r_counter) + delta_t;
AA_array(r_counter+1) = AA_array(r_counter) + chosen_delta_AA;

end

%%% Get the results into shape for plotting
% Now, get the resulting arrays into good shape
last_time = tt_array(r_counter+1);
last_step = r_counter;

if last_time > max_time
% The maximal time was reached before the maximal number of steps was
% reached. This means we have to throw away all elements of the
% initially allocated arrays, because they have not been assigned
tt_array = tt_array(1:last_step);
tt_array(end) = max_time; % Otherwise the last time will be too high
AA_array = AA_array(1:last_step);
else
disp('Not enough reaction steps allowed to reach maximal time.')
end

%%% Plot the results
plot(tt_array,AA_array,'k-')
hold on
plot(tt_array,NN-AA_array,'r--')
hold off

xlabel('Time')
ylabel('Particle count')

legend('N_A','N_B')

set(gca,'XLim',[0,tt_array(end)]) % Make axis range fit simulation

## 2 thoughts on “Simple Stochastic Simulation in MatLab”

1. An explanation of stochastic in layman’s (terms and level of knowledge) would be really well appreciated. Anyway can’t use the code without co-relating to the original understanding of the stochastic.

1. Thanks for your comment Rahul.

When I think of a stochastic process, the first property is that every time you repeat it, the result can not be predicted. You don’t know the number it will produce next.

The second property, which makes it interesting to study, is that over many repetitions, you can see patterns in the randomness.

The simplest example maybe would be tossing a coin. You never know what’s next, heads or tails. But, you do it 100 times, and you’ll see a regularity…

Now, for the example of a Gillespie process, it’s already a little more complicated. It’s a rather specific tool. You have the current state of a given system. You know which rates exist to go to different other states. Based on these rates you draw, stochastically (random, but adherent to distributions fitting the rates), the next change in the systems state, and – importantly – how long you have to wait for that change.

In this aspect it differs from many other models of stochasticity (Brownian motion/Wiener process etc.). For those you would typically approximate the process at equally spaced time points and draw a randomly distributed change from one to the next time point.

Hope that helps, good luck choosing the right type of stochastic process. There’s different kinds of random 😉