function s_arr=f_StateSequenceGillespie(Q,dt,p0,Nsteps,Nr)

% Using Gillespie's algorithm, given Markovian model "Q" 
% (single channel or composite), computes for
% "Nsteps" successive steps of duration "dt" and for "Nr" sweeps
% the state of the model at each time step and for each sweep,
% given initial probability distribution p0

% Input
% Matrix Q of the model
% dt
% Initial probability p0 as column vector (if empty [], will be set to steady state)
% Number Nsteps of time steps (output will include the first and the last value)
% Number of independent realizations Nr (sweeps)

% Output:
% s_arr: Nr-row x Nsteps+1 matrix of state indices

% Equilibrium (steady state distribution), if p0 not given explicitly
if isempty(p0)
    p0=findSS(Q);
end

% Allocate array for the result
s_arr=zeros(Nr,Nsteps+1); % includes initial and final state

for i=1:Nr

    % Initialize Gillespie's algorithm
    s=find(rand()<cumsum(p0),1); 
    stepctr=0; % step counter
    tnext=0; % Time
    
    while stepctr<Nsteps+1
       
        % Calculate the time until the next transition    
        s_outgoing=-Q(s,s); % Sum of the outgoing rates, obtained from the corresponding element of Q on the diagonal
        time=-log(1-rand())/s_outgoing;
        tnext=tnext+time; % time at which the next transition occurs

        % Calculate to which state this transition occurs    
        r=Q(:,s); % Fetch the outgoing rates
        r(s)=[]; % Remove the negative element (on the diagonal)    
        cs=cumsum(r);
        csum=cs(end);
        ns=find(rand()<cumsum(r)/csum,1);
        if ns>=s
           ns=ns+1; 
        end
        % ns is now the state to which the transition occurs
        
        % Go forward in time and insert the corresponding state in s_arr(i,:)
        % until the tramnsition occurs        
        nextstepctr=floor(tnext/dt);
        nextstepctr=min([nextstepctr Nsteps+1]);
        s_arr(i,stepctr+1:nextstepctr)=s;
                
        % update the state and the step counter
        s=ns;
        stepctr=nextstepctr;
        
    end % while

end % for
