% Fofana & Hurdford, Mechanistic movement models to understand epidemic spread,Philosophical Transactions B. DOI: 10.1098/rstb %Define spatial and temporal vectors for sim=1:30 Li=0; Lf=100; delta_x=0.1; delta_y=0.1; Lx=Li:delta_x:Lf; Ly=Li:delta_y:Lf; x_step=2; radius=1; ti=0; tf=500; tstep=1; T=ti:tstep:tf; N=1000; p_r=0.3; p_l=0.7; t=1; tupdate=1; % Set initial distribution bwtraj_x=zeros(length(T),N); bwtraj_y=zeros(length(T),N); bwpos_x=randsample(Lx,N,'true'); bwpos_y=randsample(Ly,N,'true'); bwtraj_x(t,:)=bwpos_x; bwtraj_y(t,:)=bwpos_y; % BRWs process for i=1:tf for j=1:N p=rand(1,1); angle=2*pi*rand(1,1); if pp_r bwpos_x(j)=bwpos_x(j)+cos(angle)*x_step; bwpos_y(j)=bwpos_y(j)+sin(angle)*x_step; else bwpos_x(j)=bwpos_x(j); bwpos_y(j)=bwpos_y(j); end % fix periodic boundary condition if bwpos_x(j)
  • Lf bwpos_x(j)=Lf; end if bwpos_y(j)
  • Lf bwpos_y(j)=Lf; end bwtraj_x(t+1,j)=bwpos_x(j); bwtraj_y(t+1,j)=bwpos_y(j); end t=t+tupdate; end if any(bwtraj_xLf) display('true') else display('false') end if any(bwtraj_yLf) display('true') else display('false') end bwdistances=zeros(N,N); bwcontact=zeros(N,N); bwcontacts=zeros(length(T),N); bwstatuts=zeros(length(T)+1,N); % Set initially infected individuals bwstatuts(1,1000)=1; for time=1:length(T) for focalind=1:N bwSIcont=0; loc_x1=bwtraj_x(time,focalind); loc_y1=bwtraj_y(time,focalind); for ind=1:N loc_x2=bwtraj_x(time,ind); loc_y2=bwtraj_y(time,ind); dist=sqrt((loc_x1-loc_x2)^2+(loc_y1-loc_y2)^2); bwdistances(focalind,ind)=dist; % check contact occurrence if dist<=delta_x bwcontact(focalind,ind)=1; else bwcontact(focalind,ind)=0; end % Changes the statut of infectious and newly infected infcont=bwstatuts(time,focalind)+bwstatuts(time,ind); state=bwstatuts(time,focalind); if infcont==1 && dist<=radius bwSIcont=bwSIcont+0.5; bwstatuts(time+1,focalind)=1; % generate new infections % update the state of newly infected and infectious % individuals. elseif state==1; bwstatuts(time+1,focalind)=1; end end bwSIconts(focalind)=bwSIcont; count=sum(bwcontact(focalind,:)); bwcontacts(time,focalind)=count-1; if bwcontacts(time,focalind)==-1 bwcontacts(time,focalind)=0; end end bwSIcontacts(sim,time)=sum(bwSIconts); bwsusceptible(sim,time)=histc(bwstatuts(time,:),0); bwinfected(sim,time)=histc(bwstatuts(time,:),1); end end