% an agent based simulation of SIS epidemic with diffusion in one dimension 
clear
% set some parameters
L = 100; % length of 1D domain
dt = 0.01; %the time step  

Nstep = 4000;  % the number of time steps to take
D = 1;  % the diffusion coefficient for the particles
% all species have the same diffusion coefficient
scD = sqrt(2*D*dt);  % variance for diffusion
%there are three possible reactions
rxn = 3;  % reaction (transmission) rate I+S -> 2I
% set rxn = 0 to see pure diffusion
Rxn = 1; %size of the reaction neighborhood
rmv = 0.5; % recovery rate I -> S
%rmv = 0; % for no recovery
% set rvm = 0 to see Fisher's equation in an agent based model
 

figon = 1;  %plot particle positions
% figon = 0;  % do not plot particle positions - simulation is MUCH faster   
%assign colors to the different particle types
col(1) = 'g';
col(2) = 'r';
col(3) = 'b';
 
    clear X
    clear S
   S0 = 1.1;  %This is the density of the initial susceptible population
 
% this means that there are a total of p*L particles 
N = fix(S0*L); % total number of particles
 
% randomly assign initial positions
x = rand(N,1)*L;
t = 0;
%there are two possible states
% s=1; susceptible
% s=2 ; infected
 
s = ones(N,1); %initially all are susceptible
% or all are initially infected
%s = 2*ones(N,1);

% except one infected at the center of the domain
x(N) = L/2;
s(N) = 2;
X = x; % save the positions
T (1) = t;
S = s; %Save the states;
for n = 1:Nstep

%there are two steps to take
%1) move the particles by diffusion

x = x + scD*randn(N,1);
x = abs(x); % reflecting boundaries at x=0;
x = L - abs(x-L); % reflecting boundaries at x=L;

% calculate reaction probabilities

for j = 1:N
% reaction probability for susceptibles:    
    if(s(j)==1)
        % calculate distances between susceptible and infected particles
        dist = zeros(1,N);
    for k = 1:N
        if (s(k) ==2)
        dist(k) =  (x(j)-x(k)).^2;
        end
    end
    DltR = sum((dist < Rxn^2)&(dist>0)); % find number of particles within distance Rxn
    rxp(j) = rxn*DltR/(2*Rxn);
    end
    if (s(j)==2)
        rxp(j) = rmv;% the recovery rate
    end
    
end
 
rxpp = find(rxp*dt>rand(1,N)); %These are the particles that change state
s(rxpp) = s(rxpp) + 1;  % increment the state by one
rxpp = find(s==3);  % if state was set to 3, set it back to one
s(rxpp) = 1;
% one cycle of the process is complete
t = t + dt;
T(n+1) = t;
X=[X,x];
S = [S,s];

end

figure(1)
for k = 1:N
    tp1 = find(S(k,:)==1);
    tp2 = find(S(k,:) == 2);
    tp2m= min(tp2);
 plot(T(tp1),X(k,tp1), strcat(col(1) ));
 plot(T(tp2),X(k,tp2), strcat(col(2) ),T(tp2m),X(k,tp2m),'r*');
 hold on
end
axis([0 T(end) 0 L])
xlabel('time','fontsize',20)
ylabel('x','fontsize',20)
hold off

nI = sum(S==2);


nIx = find(nI<N);
sl = sum(T(nIx).*nI(nIx))/sum(T(nIx).^2);
figure(2)
plot(T,nI/N,T(nIx),sl*T(nIx)/N,'--','linewidth',2)
xlabel('time','fontsize',20)
ylabel('% of Infected','fontsize',20)
axis([0 T(end) 0 1])

 

