% Gillespie algorithm for a simple 2 state model

% initialize the states
clear
N = 2;  % number of channels
M = 4;  % number of independent subunits per channel

alpha = 1.5; %open rate
beta = 0.5; %close rate


K = 1000; % number of time steps
T = zeros(N,1); % store the transition times
P = zeros(N,1); % store the states
p = P;
delt = [];

for j = 0:M
    rate = beta*j + (M-j)*alpha;% total rate for each state
    rates(j+1) =rate;
    ptest(j+1) = beta*j/rate; % ratio of rates
end

for j = 2:K
    % pick a random number
    k = rand(N,1);
    % use this random number to determine the next transition time
    rate = rates(p(:)+1)
    delt(:,j)=-log(k)./rate';
    T(:,j) = delt(:,j)+T(:,j-1);  %this sets the next time
   
    test = ptest(p(:)+1)'
    % pick another random number
    k = rand(N,1);
    %use this random number to determine which transition to make
    stepsize = -1+2*(k>test)
    p = p+stepsize;
    P(:,j) = p;
   
end

% for a particular channel
j = 2;
figure(1)
plot(T',P(j,:)' )

figure(2)
hist(delt(j,:),20)