
%%%%%%%%%%%%%%% Function SSA Stochastic Simulation Algorithm %%%%%%%%%%%%%%%

function ...
[X,T] = SSA(x, prop, nu, t, Kcat, Kc, CRP, Barac, Karac, Aarac, ...       
    Parac, Kc1, arabinose, Bca, Kca1, Kca2, Aca, Pca, Acamp, CAMPi, nc)

% Main loop
i=1; 
while t<400;
% 400 minutes = interesting time given by ODE modeling
    [r]=rand(1,2);   
    % Two random numbers for the SSA
    
    prop0 = sum(prop);

    % Chose tau randomly 
    t0 = (1/prop0)*log(1/r(1));

    % chose randomly the next event given the probability of each event
    j= min(find(cumsum(prop)> prop0*r(2)));

    % Update
    t = t+t0;                                                           
    T(i) = t;                                                               

    % Application of the matrix nu
    x(:) = x(:) + nu(:,j);
    
    X(i) = x(1);            
    % Concentrations of CAMP and complexes
        
        % Variable for cAMP
        CAMP = (Kcat*x(1)+0.6*CAMPi)/Acamp;                                          
        
        % Variable equivalent to the production of (CRP-cAMP)
        C1 = CRP*CAMP / (Kc + CAMP);                                            
        
        % Variable equivalent to the production of arac
        C2 = x(2) / (1 + Kc1 / arabinose);                                     

    % Propensity functions
        % 1 >> CA
        % 2 >> Arac
        
        % (CRP-cAMP) >> Arac
        prop(1) = Barac*C1^nc / (Karac^nc + C1^nc);
        
        % inhibition of arac
        prop(2) = Aarac*x(2);        
        
        % Basal production of arac
        prop(3) = Parac;
        
        % AND gate >> production of CA
        prop(4) = Bca*(C1^nc/(C1^nc+Kca1^nc))*(C2^nc / (Kca2^(nc) + C2^(nc)));   
        
        % Inhibition of CA
        prop(5) = Aca*x(1);   
        
        % Basal production of CA
        prop(6) = Pca;                                                       


    i = i+1 ;                                                              
     

end


end
    