% write an IBM for animal disease interaction %parameters %3 million paired nucleotides in average bacterial genome %change 3% for a speciation or 90,000 %the average mutation rate of a wild-type E. coli strain is ?1 × 10?3 per genome per generation - Lee et al %11,688 generations per year clear % for zzx=1:3 % if zzx==1 % disnum=0.8; % elseif zzx==2 % disnum=0.9; % else % disnum=0.95; % end chprob2 = 100; disnum=0.9; timez = 20000; spdis=0.9; spdis2=0.9; szzz=500; sza=szzz; szb=szzz; sp=zeros(sza,szb); sp2=zeros(sza,szb); dis=zeros(sza,szb);%disease pathdiff=zeros(sza,szb); ch=zeros(sza,szb);%change in genome pd1=0; tota=0; totb=0; totc=0; totd=0; %%create map of species in 5% of space for i = 50:sza-50 for j= 50:sza-50 if rand(1)>spdis % assume 5% of species randomly distributed sp(i,j)=1; end end end for i = 50:sza-50 for j= 50:sza-50 if rand(1)>spdis2 % assume 10% of species randomly distributed sp2(i,j)=1; end end end for i = 50:sza-50 for j= 50:sza-50 if rand(1)>disnum && sp(i,j)==1% assume 10% of animals with pathogen dis(i,j)=1; end if rand(1)>disnum && sp2(i,j)==1% assume 5% of species randomly distributed dis(i,j)=1; end end end %assume each time step is 10 years for z = 1:timez % imagesc(sp) % pause(0.01) for i = 5:sza-5 for j= 5:szb-5 %at each time step add one genetic counter if dis(i,j)==1 ch(i,j)=ch(i,j)+1; end %create a matrix of 9 by 9 around each species of interest if sp(i,j)==1 spt = sp(i-4:i+4,j-4:j+4); sp2t = sp2(i-4:i+4,j-4:j+4); dist = dis(i-4:i+4,j-4:j+4); cht = ch(i-4:i+4,j-4:j+4); if sum(sp2t(:))>0 && dis(i,j)==1 %if there is a gap within the dispersal distance and the tree is still alive zzz=1; for ii=1:9 for jj=1:9 if ii==5 && jj==5 else if dist(ii,jj)==1 && cht(ii,jj)>0%.98 pd1(zzz) = abs(ch(i-5+ii,j-5+jj)-ch(i,j)); ch(i-5+ii,j-5+jj)=0; % tota=tota+1; zzz=zzz+1; end end end end if nansum(pd1)>0 pathdiff(i,j) =mean(pd1); pd1=0; end end %move animal randomly by one cell zs1 = randi(2)-1; zs2 = randi(2)-1; if rand(1)>0.51 zs1 = -zs1; end if rand(1)>0.51 zs2 = -zs2; end %first move animal without disease if dis(i,j)==0 if sp(i+zs1,j+zs2)==1 || sp2(i+zs1,j+zs2)==1 sp(i,j)=1; else sp(i,j)=0; sp(i+zs1,j+zs2)=1; end else %then move animal with disease if sp(i+zs1,j+zs2)==1 || sp2(i+zs1,j+zs2)==1 else sp(i,j)=0; sp(i+zs1,j+zs2)=1; if dis(i,j)==1 dis(i,j)=0; dis(i+zs1,j+zs2)=1; elseif dis(i,j)==2 dis(i,j)=0; dis(i+zs1,j+zs2)=2; elseif dis(i,j)==0 end ch(i+zs1,j+zs2)=ch(i,j); ch(i,j)=0; end end end %move the next species zs1 = randi(2)-1; zs2 = randi(2)-1; if rand(1)>0.51 zs1 = -zs1; end if rand(1)>0.51 zs2 = -zs2; end if sp2(i,j)==1 %first move animal without disease if dis(i,j)==0 if sp2(i+zs1,j+zs2)==1 || sp(i+zs1,j+zs2)==1 sp2(i,j)=1; else sp2(i,j)=0; sp2(i+zs1,j+zs2)=1; end else %then move animal with disease if sp2(i+zs1,j+zs2)==1 || sp(i+zs1,j+zs2)==1 else sp2(i,j)=0; sp2(i+zs1,j+zs2)=1; if dis(i,j)==1 dis(i,j)=0; dis(i+zs1,j+zs2)=1; elseif dis(i,j)==2 dis(i,j)=0; dis(i+zs1,j+zs2)=2; elseif dis(i,j)==0 end ch(i+zs1,j+zs2)=ch(i,j); ch(i,j)=0; end end end end end chx(z)= length(find(dis==2)); ch1x(z)= length(find(dis==1)); spx(z)= length(find(sp==1)); sp2x(z)= length(find(sp2==1)); chzz = ch(100:sza-100,100:szb-100); chzz2 = chzz(:); mr1(z) = sum(chzz2)/sum(dis(:)); mrxx1(z) = max(chzz2); mrx1(z) = length(find(chzz2>chprob2))/sum(dis(:)); pzz = pathdiff(100:sza-100,100:szb-100); pzz2 = pzz(:); pzz3 = pzz2(find(pzz2>0)); pdx1(z) = mean(pzz3); pdx1med(z) = median(pzz3); pdx1max(z) = max(pzz3); if timez==timez(end) pathdiff1=pathdiff; else pathdiff=zeros(sza,szb); end end disz=dis; h1=ch(:); length(find(h1>1)) %plot(ch1x,'k') clear ch sp=zeros(sza,szb); sp2=zeros(sza,szb); ch=zeros(sza,szb);%change in genome dis=zeros(sza,szb);%disease dis2=zeros(sza,szb);%disease ch=zeros(sza,szb);%change in genome pathdiff=zeros(sza,szb); %%create map of species in 5% of space for i = 50:sza-50 for j= 50:szb-50 if rand(1)>spdis % assume 5% of species randomly distributed sp(i,j)=1; end end end for i = 50:sza-50 for j= 50:szb-50 if rand(1)>spdis2 % assume 10% of species randomly distributed sp2(i,j)=1; end end end for i = 50:sza-50 for j= 50:szb-50 if rand(1)>disnum && sp(i,j)==1% assume 10% of animals with pathogen dis(i,j)=1; end if rand(1)>disnum && sp2(i,j)==1% assume 5% of species randomly distributed dis(i,j)=1; end end end %for reduced home range for z = 1:timez for i = 5:sza-5 for j= 5:szb-5 if dis(i,j)==1 ch(i,j)=ch(i,j)+1; end %create a matrix of 9 by 9 around each species of interest if sp(i,j)==1 spt = sp(i-1:i+1,j-1:j+1); sp2t = sp2(i-1:i+1,j-1:j+1); dist = dis(i-1:i+1,j-1:j+1); cht = ch(i-1:i+1,j-1:j+1); if sum(sp2t(:))>0 && dis(i,j)==1 %if there is a gap within the dispersal distance and the tree is still alive zzz=1; for ii=1:3 for jj=1:3 if ii==2 && jj==2 else if dist(ii,jj)==1 && cht(ii,jj)>0 %.98 pd1(zzz) = abs(ch(i-2+ii,j-2+jj)-ch(i,j)); ch(i-2+ii,j-2+jj)=0; % totb=totb+1; zzz=zzz+1; end end end end if nansum(pd1)>0 pathdiff(i,j) =mean(pd1); pd1=0; end end %move animal randomly by one cell zs1 = randi(2)-1; zs2 = randi(2)-1; if rand(1)>0.51 zs1 = -zs1; end if rand(1)>0.51 zs2 = -zs2; end %first move animal without disease if dis(i,j)==0 if sp(i+zs1,j+zs2)==1 || sp2(i+zs1,j+zs2)==1 sp(i,j)=1; else sp(i,j)=0; sp(i+zs1,j+zs2)=1; end else %then move animal with disease if sp(i+zs1,j+zs2)==1 || sp2(i+zs1,j+zs2)==1 else sp(i,j)=0; sp(i+zs1,j+zs2)=1; if dis(i,j)==1 dis(i,j)=0; dis(i+zs1,j+zs2)=1; elseif dis(i,j)==2 dis(i,j)=0; dis(i+zs1,j+zs2)=2; elseif dis(i,j)==0 end ch(i+zs1,j+zs2)=ch(i,j); ch(i,j)=0; end end end %move the next species zs1 = randi(2)-1; zs2 = randi(2)-1; if rand(1)>0.51 zs1 = -zs1; end if rand(1)>0.51 zs2 = -zs2; end if sp2(i,j)==1 %first move animal without disease if dis(i,j)==0 if sp2(i+zs1,j+zs2)==1 || sp(i+zs1,j+zs2)==1 sp2(i,j)=1; else sp2(i,j)=0; sp2(i+zs1,j+zs2)=1; end else %then move animal with disease if sp2(i+zs1,j+zs2)==1 || sp(i+zs1,j+zs2)==1 else sp2(i,j)=0; sp2(i+zs1,j+zs2)=1; if dis(i,j)==1 dis(i,j)=0; dis(i+zs1,j+zs2)=1; elseif dis(i,j)==2 dis(i,j)=0; dis(i+zs1,j+zs2)=2; elseif dis(i,j)==0 end ch(i+zs1,j+zs2)=ch(i,j); ch(i,j)=0; end end end end end ch2x(z)= length(find(dis==2)); ch1x(z)= length(find(dis==1)); sp2x(z)= length(find(sp==1)); sp22x(z)= length(find(sp2==1)); chzz = ch(100:sza-100,100:szb-100); chzz2 = chzz(:); mr2(z) = sum(chzz2)/sum(dis(:)); mrxx2(z) = max(chzz2); mrx2(z) = length(find(chzz2>chprob2))/sum(dis(:)); pzz = pathdiff(100:sza-100,100:szb-100); pzz2 = pzz(:); pzz3 = pzz2(find(pzz2>0)); pdx2(z) = mean(pzz3); pdx2med(z) = median(pzz3); if nansum(pzz3)~=0 pdx2max(z) = max(pzz3); end if timez==timez(end) pathdiff2=pathdiff; else pathdiff=zeros(sza,szb); end end h1=ch(:); length(find(h1>1)); dis2z=dis; subplot(311) plot(pdx1,'r') hold on plot(pdx2) fame(zzx) = pdx2(end)/pdx1(end); famed(zzx) = mrx2(end)/mrx1(end); famd(zzx) = mrxx2(end)/mrxx1(end); zzzz=[fame; famed; famd]; subplot(312) plot(mrx1,'r') hold on plot(mrx2) subplot(313) plot(mrxx1,'r') hold on plot(mrxx2) %end %