% 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_l=0.5; p_r=0.5; t=1; tupdate=1; % Set initial distribution rwtraj_x=zeros(length(T),N); rwtraj_y=zeros(length(T),N); rwpos_x=randsample(Lx,N,'true'); rwpos_y=randsample(Ly,N,'true'); rwtraj_x(t,:)=rwpos_x; rwtraj_y(t,:)=rwpos_y; % URW process for i=1:tf for j=1:N p=rand(1,1); angle=2*pi*rand(1,1); if pp_l rwpos_x(j)=rwpos_x(j)+cos(angle)*x_step; rwpos_y(j)=rwpos_y(j)+sin(angle)*x_step; else rwpos_x(j)=rwpos_x(j); rwpos_y(j)=rwpos_y(j); end % fix periodic boundary condition if rwpos_x(j)
  • Lf rwpos_x(j)=Lf; end if rwpos_y(j)
  • Lf rwpos_y(j)=Lf; end rwtraj_x(t+1,j)=rwpos_x(j); rwtraj_y(t+1,j)=rwpos_y(j); end t=t+tupdate; end if any(rwtraj_xLf) display('true') else display('false') end if any(rwtraj_yLf) display('true') else display('false') end % the contact and infection process. rwdistances=zeros(N,N); rwcontact=zeros(N,N); rwcontacts=zeros(length(T),N); rwstatuts=zeros(length(T)+1,N); % Set one individual infectious at the initial time (0 is uninfected and 1 is infected) rwstatuts(1,1000)=1; for time=1:length(T) for focalind=1:N rwSIcont=0; loc_x1=rwtraj_x(time,focalind); loc_y1=rwtraj_y(time,focalind); for ind=1:N loc_x2=rwtraj_x(time,ind); loc_y2=rwtraj_y(time,ind); dist=sqrt((loc_x1-loc_x2)^2+(loc_y1-loc_y2)^2); rwdistances(focalind,ind)=dist; % check rwcontact occurrence if dist<=radius rwcontact(focalind,ind)=1; else rwcontact(focalind,ind)=0; end % check infection occurrence % Changes the statut of infectious and newly infected infcont=rwstatuts(time,focalind)+rwstatuts(time,ind); state=rwstatuts(time,focalind); if infcont==1 && dist<=radius rwSIcont=rwSIcont+0.5; rwstatuts(time+1,focalind)=1; % generate new infections and update the state of newly % infected individuals. % Then update the state of infected individuals elseif state==1; rwstatuts(time+1,focalind)=1; end end rwSIconts(focalind)=rwSIcont; count=sum(rwcontact(focalind,:)); rwcontacts(time,focalind)=count-1; if rwcontacts(time,focalind)==-1 rwcontacts(time,focalind)=0; end end rwSIcontacts(sim,time)=sum(rwSIconts); rwsusceptible(sim,time)=histc(rwstatuts(time,:),0); rwinfected(sim,time)=histc(rwstatuts(time,:),1); end end