function [X, t] = ssa(c,t_fin,X0)

% function for the SSA for simulation of GRNs  written by Kevin Burrage
% (UQ) and Margherita Carletti (Urbino) - the user has to change the
% stoichiometric vectors and the propensity functions.
% There is a clever sort procedure to find the reaction number for a given
% step.
% The rate constants are in the c vector passed in from the main program


% time interval [0, t_fin]
% initial condition  X0

% state change matrix
v1=[0 0 1 0 0 0]';
v2=[1 0 0 0 0 0]'; v3=[-1 0 0 0 0 0]'; v4=[0 0 -1 0 0 0]'; v5=[-2 1 0 0 0 0]'; 
v6=[2 -1 0 0 0 0]'; v7=[0 -1 0 -1 1 0]'; v8=[0 1 0 1 -1 0]'; v9=[0 -1 0 0 -1 1]'; v10=[0 1 0 0 1 -1]';

v=[v1 v2 v3 v4 v5 v6 v7 v8 v9 v10];
rand('state',0);


t_old  = 0;   
x      = X0';
t(1)   = 0;
X(:,1) = X0;

count  = 1;        

while (t_old  < t_fin)
      a(1) = c(1)*x(5);
      a(2) = c(2)*x(3);
      a(3) = c(3)*x(1);
      a(4) = c(4)*x(3);
      a(5) = c(5)*x(1)^2;
      a(6) = c(6)*x(2);
      a(7) = c(7)*x(2)*x(4);
      a(8) = c(8)*x(5);
      a(9) = c(9)*x(2)*x(5);
      a(10)= c(10)*x(6);
      a=[a(1) a(2) a(3) a(4) a(5) a(6) a(7) a(8) a(9) a(10)]; 
      a_0 = sum(a);
      if a_0==0
          disp('all populations go to extinction...program stops here')
          break
      end;
       
      tau = 1/a_0 * log(1/rand);
      standard = rand * a_0;
      A=cumsum(a);
      A=[A standard];
      B=sort(A);
      pos=find(B==standard);
      x=x+v(:,pos)';
      t_old=t_old+tau;
      t(count) = t_old;
      X(:,count) = x';
      count = count + 1;
    
end
a
t(count) = t_old;
X(:,count) = x';
size(t)
size(X)