Gillespie's Algorithm (SSA)
This script implements Gillespie's Algorithm taking advantage of Matlab's vector operations. The for loops have been kept at a minimum. Gillespie's Algorithm is considered an exact stochastic simulation, however the downside is that it is much more computationally expensive than other approximate methods. This is why so much effort was put into vectorizing the script. In this model we have four species of partivles and eight possible reations that can occur. The method works by using a logarithmic distribution to find the time until the next reaction occurs, and then randomly picks from the reactions which occurs and then updates the species. For more information on the method please see (Gillespie, 1977). The parameters here are based off of Morshed et al., 2017 in BioSystems.
close all;clear; rng ('shuffle');
We have four species, eight reactions, an ensemble of 10000 and will let it run for 100 reactions.
n_species = 4; n_paths = 8; n_react = 1e2; n_ens = 100;
epsilon = 0.05; stiff = 1000; alpha1 = 1.26; alpha2 = 1.26; beta = 4; gamma = 4; ii1 = 0; ii2 = 0; c(1)=alpha1*23; c(2)=0.23; c(3)=0.23; c(4)=0.23; c(5)=alpha2*23; c(6)=0.23; c(7)=0.23; c(8)=0.23;
This defines the results of the various reactions on the species.
react_result = [0 0 1 0; 0 0 -1 0; 1 0 0 0; -1 0 0 0; 0 0 0 1; 0 0 0 -1; 0 1 0 0; 0 -1 0 0];
Initial numbers set to zero
X0 = [76;75;60;60]; n_mol = repmat(X0,1,n_ens); prop = zeros(n_paths,n_ens); react_result=repmat(react_result,1,1,n_ens); ti=zeros(1,n_ens); tn=ti; cntr=0; for ii=1:n_react
The propensity of each reaction.
prop(1,:) = (stiff*c(1)*(1+ii2)^beta)./(1+ii2+0.01*n_mol(2,:)).^beta; prop(2,:) = stiff*c(2)*n_mol(3,:); prop(3,:) = c(3)*n_mol(3,:); prop(4,:) = c(4)*n_mol(1,:); prop(5,:) = stiff*(c(5)./(1+((0.01*n_mol(1,:))./(1+ii1)).^gamma)); prop(6,:) = stiff*c(6)*n_mol(4,:); prop(7,:) = c(7)*n_mol(4,:); prop(8,:) = c(8)*n_mol(2,:);
Time to next reaction
deltat=-log(rand(1,n_ens))./sum(prop); tn = tn + deltat;
See valid reactions
The greater than returns a logical matrix which I then element-wise multiply to keep the matrix the same size.
valid_ind = prop>0; valid_prop = prop.*valid_ind; valid_ind_3d = repmat(valid_ind,1,1,n_species); valid_ind_3d = reshape(valid_ind_3d,[n_paths n_species n_ens]); valid_result = react_result.*valid_ind_3d;
Record the locations of zeros.
I take advantage that the default for most operations is to apply them column-wise.
prop_int = cumsum(valid_prop); % Set the invalid paths to NaN to not mess anything up. prop_int(invalid)=NaN;
Use bsxfun to normalize each row (species) of the propensity intervals by a vector that is the maximum of each column,
prop_int = bsxfun(@rdivide,prop_int,max(prop_int));
Pick the reaction
Make the random vector with one number for each ensemble member.
rand_react = rand(1,n_ens);
Find all propensity intervals that are greater than the random number for each ensemble (column).
react_ind = bsxfun(@gt,prop_int,rand_react);
Find the index of the first entry in each column that is nonzero.
[~,ind] = max(react_ind~=0,,1);
Since this gives the respective row index that the first nonzero entry appears we need to multiply by the number of species times the ensemble number to get the position in linear indexing. Had to make this more complicated for the general case. This is the starting index for each of the corresponding rows in linear indexing.
start_ind = (0:n_ens-1).*n_species.*n_paths + (1 + n_species.*(ind-1));
Since we want the entire 'row' or all the species
species = 0:n_species-1;
Add replicated versions of the two vectors to endup with all the relavent indices.
indicies = bsxfun(@plus,repmat(start_ind,n_species,1),repmat(species(:),1,n_ens));
Add the valid result for to each species for the current number of mols.
n_mol = n_mol+valid_result(indicies);
figure(1) clf subplot(2,2,1) hist(n_mol(1,:)) subplot(2,2,2) hist(n_mol(2,:)) subplot(2,2,3) hist(n_mol(3,:)) subplot(2,2,4) hist(n_mol(4,:))