%Stochastic SIR frequency-dependent model

tmax=1; %length of time (in years) to simulate

R0=10;
gamma=365/10;
mu=1/70;
beta=R0*(gamma+mu);

Ninit=10000;

Sinit=1500; %initial susceptible population
Iinit=1; %initial infected population

% Code to initialize at the deterministic equilibrium
% Sinit=floor(Ninit/R0); %initial susceptible population
% Iinit=floor(Ninit*(R0-1)*mu/beta); %initial infected population

Rinit=Ninit-Iinit-Sinit; %initial recovered/immune population

nsims=10; %number of simulations to run

for i=1:nsims   %simulation loop
    clear time suscep infect recover pop
    j=1;
    t=0;
    time(j)=t;
    S=Sinit; 
    I=Iinit; 
    R=Rinit; 
    N=S+I+R;
    suscep(j)=S;
    infect(j)=I;
    recover(j)=R;    
    pop(j)=N;
    % Events
    a1=mu*N;
    a2=mu*S;
    a3=beta*S*I/N;
    a4=mu*I;
    a5=gamma*I;
    a6=mu*R;
    while (t<tmax)
        A=a1+a2+a3+a4+a5+a6; % sum of all events 
        y1=rand;             % generates random number between zero and one
        dt=-log(y1)/A;       % use random number to estimate time to next event 
        y2=rand;             % use random number to decide which event happens
        if y2<(a1/A)                      
            S=S+1; 
            N=N+1;
            a1=mu*N;
            a2=mu*S;
        elseif y2<((a1+a2)/A)
            S=S-1;
            N=N-1;
            a1=mu*N;
            a2=mu*S;
        elseif y2<((a1+a2+a3)/A)
            S=S-1;
            I=I+1;
            a2=mu*S;
            a4=mu*I;
            a5=gamma*I;
        elseif y2<((a1+a2+a3+a4)/A)
            I=I-1;
            N=N-1;
            a4=mu*I;
            a5=gamma*I;
        elseif y2<((a1+a2+a3+a4+a5)/A) 
            I=I-1;
            R=R+1;
            a4=mu*I;
            a5=gamma*I;
            a6=mu*R;
        else
            R=R-1;
            N=N-1;
            a6=mu*R;
        end; 
        a3=beta*S*I/N;            
        t=t+dt;
        j=j+1;
        time(j)=t; 
        suscep(j)=S;
        infect(j)=I;
        recover(j)=R;
        pop(j)=N;  
    end; %while loop
%     subplot(2,1,1);
%     plot (time,[pop;suscep;recover])
%     ylabel('Number');
%     legend('N','S','R')
%     subplot(2,1,2);
    plot(time,infect)
    ylabel('Number of Infectives');
    xlabel('Time (years)');
    hold on
%    pause
end; %simulation loop

    



