% A Two-Patch SIS program for MBI-NIMBioS-CAMBAM PROGRAM
% Continuous-time Markov chain model
% nsamp sample paths are calculated 
% s1,s2,i1,i2 are population sizes in two patches
% gam1,gam2,di12,di21,ds12,ds21,bet1,bet2 are constants
% s10,s20,i10,i20 are the initial population sizes
clear all
tmax=25.; %***
nsamp=1000; %***
%Initial conditions
s10=25; s20=20; %***
i10=30; i20=25; %*** 
%Parameter values
gam1=.05; gam2=.4; %*** 
di12=.3; di21=.05; %*** 
ds12=0; ds21=.25; %*** 
bet1=.5; bet2=.01; %***
%Patch reproduction numbers
r01=bet1/gam1;
r02=bet2/gam2;
j=1;
i1(j)=i10; i2(j)=i20;
s1(j)=s10; s2(j)=s20;
tot1=s1(j)+i1(j); tot2=s2(j)+i2(j);
tt(j)=0;
%One sample path
while tt(j)<tmax & tot1*tot2>0;
    ra1=rand;
    ra2=rand;
    tot1=s1(j)+i1(j);
    tot2=s2(j)+i2(j);
    ev1=bet1*s1(j)*i1(j)/tot1;
    ev2=bet2*s2(j)*i2(j)/tot2;
    ev3=gam1*i1(j);
    ev4=gam2*i2(j);
    ev5=ds12*s1(j);
    ev6=ds21*s2(j);
    ev7=di12*i1(j);
    ev8=di21*i2(j);
    totev=ev1+ev2+ev3+ev4+ev5+ev6+ev7+ev8;
    tt(j+1)=tt(j)-log(ra1)/totev;
    if ra2<ev1/totev
        s1(j)=s1(j)-1;
        i1(j)=i1(j)+1;
    elseif ra2>=ev1/totev & ra2<(ev1+ev2)/totev
        s2(j)=s2(j)-1;
        i2(j)=i2(j)+1;
    elseif ra2>=(ev1+ev2)/totev & ra2<(ev1+ev2+ev3)/totev
        s1(j)=s1(j)+1;
        i1(j)=i1(j)-1;
    elseif ra2>=(ev1+ev2+ev3)/totev & ra2<(ev1+ev2+ev3+ev4)/totev
        s2(j)=s2(j)+1;
        i2(j)=i2(j)-1;
    elseif ra2>=(ev1+ev2+ev3+ev4)/totev & ra2 <(ev1+ev2+ev3+ev4+ev5)/totev
        s1(j)=s1(j)-1;
        s2(j)=s2(j)+1;
    elseif ra2>=(ev1+ev2+ev3+ev4+ev5)/totev & ra2 <(ev1+ev2+ev3+ev4+ev5+ev6)/totev
        s1(j)=s1(j)+1;
        s2(j)=s2(j)-1; 
    elseif ra2>=(ev1+ev2+ev3+ev4+ev5+ev6)/totev & ra2 <(ev1+ev2+ev3+ev4+ev5+ev6+ev7)/totev
        i1(j)=i1(j)-1;
        i2(j)=i2(j)+1;
    elseif ra2>(ev1+ev2+ev3+ev4+ev5+ev6+ev7)/totev & ra2<=1
        i1(j)=i1(j)+1;
        i2(j)=i2(j)-1;
    end
    s1(j+1)=s1(j); 
    s2(j+1)=s2(j);
    i1(j+1)=i1(j);
    i2(j+1)=i2(j);
    j=j+1;
end
subplot(2,2,1);
plot(tt,s1,'k-',tt,i1,'r--','linewidth',2);
xlabel('Time t');
ylabel('S_1 and I_1')
legend('S_1', 'I_1')
title('CTMC')
axis([0,tmax,0,60]);
hold off
subplot(2,2,2);
plot(tt,s2,'k-',tt,i2,'r--','linewidth',2);
xlabel('Time t');
ylabel('S_2 and I_2')
title('CTMC');
axis([0,tmax,0,60]);
legend('S_2', 'I_2');
hold off
clear s1 s2 i1 i2 tt
%Many sample paths
for k=1:nsamp
    j=1;
    i1(1)=i10;
    i2(1)=i20;
    s1(1)=s10;
    s2(1)=s20;
    tot1=s1(1)+i1(1);
    tot2=s2(1)+i2(1);
    tt(1)=0;
while tt(j)<tmax & tot1*tot2>0;
    ra1=rand;
    ra2=rand;
    tot1=s1(1)+i1(1);
    tot2=s2(1)+i2(1);
    ev1=bet1*s1(1)*i1(1)/tot1;
    ev2=bet2*s2(1)*i2(1)/tot2;
    ev3=gam1*i1(1);
    ev4=gam2*i2(1);
    ev5=ds12*s1(1);
    ev6=ds21*s2(1);
    ev7=di12*i1(1);
    ev8=di21*i2(1);
    totev=ev1+ev2+ev3+ev4+ev5+ev6+ev7+ev8;
    tt(j+1)=tt(j)-log(ra1)/totev;
    if ra2<=ev1/totev
        s1(1)=s1(1)-1;
        i1(1)=i1(1)+1;
    elseif ra2>ev1/totev & ra2<=(ev1+ev2)/totev
        s2(1)=s2(1)-1;
        i2(1)=i2(1)+1;
    elseif ra2>(ev1+ev2)/totev & ra2<=(ev1+ev2+ev3)/totev
        s1(1)=s1(1)+1;
        i1(1)=i1(1)-1;
    elseif ra2>(ev1+ev2+ev3)/totev & ra2<=(ev1+ev2+ev3+ev4)/totev
        s2(1)=s2(1)+1;
        i2(1)=i2(1)-1;
    elseif ra2>(ev1+ev2+ev3+ev4)/totev & ra2<=(ev1+ev2+ev3+ev4+ev5)/totev
        s1(1)=s1(1)-1;
        s2(1)=s2(1)+1;
    elseif ra2>(ev1+ev2+ev3+ev4+ev5)/totev & ra2<=(ev1+ev2+ev3+ev4+ev5+ev6)/totev
        s1(1)=s1(1)+1;
        s2(1)=s2(1)-1; 
    elseif ra2>(ev1+ev2+ev3+ev4+ev5+ev6)/totev & ra2<=(ev1+ev2+ev3+ev4+ev5+ev6+ev7)/totev
        i1(1)=i1(1)-1;
        i2(1)=i2(1)+1;
    elseif ra2>(ev1+ev2+ev3+ev4+ev5+ev6+ev7)/totev & ra2<=1
        i1(1)=i1(1)+1;
        i2(1)=i2(1)-1;
    end
    
    s1(2)=s1(1); 
    s2(2)=s2(1);
    i1(2)=i1(1);
    i2(2)=i2(1);
    s1(1)=s1(2);
    s2(1)=s2(2);
    i1(1)=i1(2);
    i2(1)=i2(2);
    j=j+1;

end
    s1e(k)=s1(1);
    s2e(k)=s2(1);
    i1e(k)=i1(1);
    i2e(k)=i2(1);
end
%Probability histograms for I1 and I2
subplot(2,2,3)
hist(i1e,[0:1:100]);
axis([-1,60,0,100]);
xlabel('n');
% nsamp=10^3
ylabel('Prob\{I_1=n\}\times 10^3');
hold off
subplot(2,2,4)
hist(i2e,[0:1:100]);
axis([-1,60,0,100]);
xlabel('n');
ylabel('Prob\{I_2=n\}\times 10^3');
prob0=sum(i1e+i2e==0)/nsamp;

%Display the means and standard deviations at time tmax
disp('CTMC  Calculational Results (Means and Standard Deviations)');
disp('     S1(mean)    I1(mean)     S2(mean)     I2(mean)') 
disp((sprintf(' %12.6f %12.6f %12.6f %12.6f', mean(s1e),  mean(i1e), mean(s2e), mean(i2e))));
disp('      S1(std)     I1(std)     S2(std)      I2(std)') 
disp((sprintf(' %12.6f %12.6f %12.6f %12.6f', std(s1e), std(i1e), std(s2e), std(i2e))));
disp('Prob Epidemic Ends    Time');
disp((sprintf('  %12.6f     % 12.6f', prob0,  tmax)));
disp('Patch 1 Reprod No.   Patch 2 Reprod No.');
disp((sprintf('%15.6f %15.6f',   r01,   r02)));
hold off