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. , with forward reaction () rate and reverse reaction () rate . Additionally a constant number of particles, $A+B=N$, was assumed, allowing the replacement , 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