gillispie_CGE.m

function stateMat = gillispie_CGE(numSteps,totalTime)
%gillespie algorithm for Constitutuve Gene Experssion (page294 of Ingalls)
% R1: null -> M    ; propensity = kr
% R2: null -> P    ; propensity = kp * nM
% R3: M    -> null ; propensity = dr * nM
% R4: P    -> null ; propensity = dp * nP
nC = 2 ; % number of chemical components
nR = 4 ; % number of chemical reactions
kr = 10; kp = 6; dr = 1; dp = 1;
% starting condition
nM = 0; nP = 0;
% state vector = [nM; nP]  (with nM = RNA and nP = Protein)
state = [nM nP];
sMatrix = [ 1 0 ; 0 1 ; -1 0 ; 0 -1];

time = 0;
step = 0;

for step = 1:1:numSteps
%while (time < totalTime)
    % calculate propensities
    nM = state(1);
    nP = state(2);
    propensitiesVec = [kr   kp*nM   dr*nM   dp*nP]
    
    % calculate totalPropensity
    totalPropensity = sum(propensitiesVec);
    
    % calculate waitTime 
    waitTime = - log( rand(1,1) ) / totalPropensity 
    
    % calculate nextReaction
    
    %prob = rand(1,1);
    %rxn = 0;
    %p1 = 0;
    %while (prob > 0.0 )
    %    rxn = rxn + 1 ;
    %    prob = prob -  propensitiesVec(rxn);
    %end
    %nextReaction = rxn;
    
    prob = rand(1,1)*totalPropensity ;
    cumulativePropensity = cumsum(propensitiesVec) ;
    rxn = 1 ;
    while ( prob > cumulativePropensity(rxn))
        rxn = rxn + 1 ;
    end
    nextReaction = rxn ;
                
    % perform nextReaction on state vector
    state = state + sMatrix(nextReaction,:);
    time = time + waitTime;
    step = step +1;
    %propensitiesVec
    stateMat(step,:) = [time, state, nextReaction];

end