## Reactivation of latent infections with migration shapes population-level disease dynamics ## last updated 07/23/2020 ## written by daniel becker ## danbeck@iu.edu ## ready workspace rm(list=ls()) graphics.off() ## load libraries library(ggplot2) library(deSolve) library(plyr) library(tidyr) ## define ode for infection dynamics at the breeding grounds mig_sili_breed=function(t,y,params){ ## state variables S=y[1] ## migrant S I=y[2] ## migrant I L=y[3] ## migrant L ## outline ODEs with( as.list(params),{ ## susceptible dS=((b0-b1*(S+I+L))*(Sbreed+(Ibreed*(1-cf))+Lbreed)*exp(-mu_b*(t-bstart))) - (mu_b*S) - (beta*S*I) ## infected dI=(beta*S*I) - (mu_b*I/(1-cs)) - (rho*I) ## latent dL=(rho*I) - (mu_b*L) ## combine dy=c(dS,dI,dL) list(dy) } ) } ## define ode for infection dynamics at the wintering grounds mig_sili_w=function(t,y,params){ ## state variables S=y[1] ## migrants S I=y[2] ## migrants I L=y[3] ## migrants L ## outline ODEs with( as.list(params),{ ## susceptible dS= - (mu_nb*S) - (beta*S*I) ## infected dI=(beta*S*I) - (mu_nb*I/(1-cs)) - (rho*I) ## latent dL=(rho*I) - (mu_nb*L) ## combine dy=c(dS,dI,dL) list(dy) } ) } ## define ode for infection dynamics during migration mig_sili_mig=function(t,y,params){ ## state variables S=y[1] ## migrants S I=y[2] ## migrants I L=y[3] ## migrants L ## outline ODEs with( as.list(params),{ ## susceptible dS= - (mu_m*S) - (beta*S*I) ## infected dI=(beta*S*I) - (mu_m*I/(1-cs)) - (rho*I) ## latent dL=(rho*I) - (mu_m*L) ## combine dy=c(dS,dI,dL) list(dy) } ) } ## divide year into timestep taustep=1000 ## takes far longer to run, only for smoothing taustep=200 ## faster but still computationally expensive ts=(1/taustep)*(1:taustep) ## within-year timesteps ## rounding factors lstep=log10(taustep) ## winter depart in March winter_depart=round(2.51/12,lstep) ## mid March ## migration as 1.5-2 months (March 1 to May 1) mig=round(1.5/12,lstep) ## 1.5 months ## arrive at breeding site b_start=round(winter_depart+mig,lstep) ## depart breeding October 1 b_depart=round(9/12,lstep) ## winter arrival winter_arrive=round(b_depart+mig,lstep) ## check all timepoints match the within-year timestep tp=c(winter_depart,b_start,b_depart,winter_arrive) tp%in%ts ## length of winter Tmw=(1-winter_arrive)+winter_depart Tmw*12 ## estimate length of the breeding season Tmb=b_depart-b_start Tmb*12 ## migration mortality (Ketterson & Nolan 1982) sigma_m=1-0.21 mu_m=-log(sigma_m)/mig ## max per capita number of offspring (Nolan et al. 2002) mF=(2*4)/2 ## estimate b0 as mF = exp(b0/t), where t = 1/t of the year being the breeding season t=1/Tmb b0=log(mF)*t ## lower bound of reproduction if 1 brood (Nolan et al. 2002) b0min=log(mF/2)*t ## K (for deriving b1) K=1000 ## survival probability during breeding (Ketterson & Nolan 1982) sigma_b=1-0.05 mu_b=-log(sigma_b)/(Tmb) ## winter survival probability (Ketterson & Nolan 1982) sigma_w=1-mean(c(0.15)) mu_nb=-log(sigma_w)/Tmw ## initial conditions N0=K S0=N0-1 I0=1 L0=0 ## cost of infection for fecundity cf=0.2 ## as in Hall et al 2014 ## cost of infection for survival (non-migration) cn=cf ## maximum years of simulation ymax=100 ## cutoff for ending the while loop ctf=0.0001 ## function to vary model parameters par_vary=function(pvar, ## proportion relapse tseason, ## transmission season rvar, ## 1/duration of infection bvar, ## transmission rate cvar, ## cost of infection for survival for migrants cnvar, ## cost of infection for survival (B and NB) cfvar, ## cost of infection for fecundity smvar, ## % survival during migration swvar, ## % survival during nonbreeding b0var, ## per capita fecundity during breeding rtime ## relapse at start or end ){ ## set pars per value rprop=pvar rho=rvar beta=bvar cs=cvar cn=cnvar cf=cfvar sigma_m=smvar sigma_w=swvar b0=b0var relapse=rtime ## convert sigma to mu mu_m=-log(sigma_m)/mig mu_nb=-log(sigma_w)/Tmw ## obtain b1 based on b0 and K b1=b0/K ## set transmission season beta_mig=ifelse(any(tseason=="mig"),beta,0) beta_breed=ifelse(any(tseason=="breed"),beta,0) beta_winter=ifelse(any(tseason=="winter"),beta,0) ## output vector giving the population sizes at each time step output=cbind(0,S0,I0,L0) ## make empty prevalence matrix pmat=matrix(nrow=taustep,ncol=1) ## fill with starting values pmat[,1]=I0/(S0+I0+L0) ## initiate while loop Y=1 ## run for 25 years, then look at the difference in each timepoint over the 2 years while(ifelse(Y<25,TRUE, ifelse(max(abs(apply(pmat[,(Y-(2-1)):Y],1,function(x) max(abs(diff(x)))))) > ctf, Y= winter_arrive){ ## wintering grounds ## relapse upon arrival after fall migration if(isTRUE(all.equal(tau,winter_arrive)) & relapse=='end'){ ## get pops at very end of migration ystart=tail(output,1)[2:4] Sm=ystart[1] Im=ystart[2] Lm=ystart[3] ## % percent of L relapse to I move=Lm*rprop ## move Lm=Lm-move Im=Im+move ## starting values ystart=c(Sm,Im,Lm) ## else just use those starting values }else{ ## starting ystart=tail(output,1)[2:4] } ## set pars pars=c(mu_nb=mu_nb,rho=rho, cs=cn,beta=beta_winter) ## set times times=c(tau,tau+((1/taustep))) # within year timestep ## solve ode answer=lsoda(ystart,times,mig_sili_w,pars) ## save output=rbind(output,c(timecounter,tail(answer,n=1)[2:4])) } else if(tau>=winter_depart & tau=b_start & tau=b_depart & tau winter_arrive,"wintering", ifelse(set$annual >= winter_depart & set$annual =b_start & set$annual0,1,0) set2$mig_trans=ifelse(beta_mig>0,1,0) set2$winter_trans=ifelse(beta_winter>0,1,0) ## make unique combination of parameters set2$upars=with(set2,paste(rprop,rho,beta, cs,sigma_m,sigma_w,cf,cf, b0,b1,relapse)) ## wide to long set4=gather(set2,group,prev,S:L) set4$group=factor(set4$group) ## reorder and rename set4$group2=revalue(set4$group,c("I"="infected", "L"="latent", "S"="susceptible")) set4$group2=factor(set4$group2,levels=c("susceptible","infected","latent")) ## get final two years final2=set2[set2$years%in%tail(unique(set2$years),2),] ## difference pm=matrix(final2$I_prop,nrow=taustep,ncol=2,byrow=F) ts=max(abs(apply(pm,1,function(x) max(abs(diff(x)))))) ## are the maxs more or less equal cycle=as.character(ifelse(ts[1]!=ts[2],'cycle','steady')) cycle=ifelse(ts < 0.001,'steady','cycle') ## add cycle to sets set2$cycle=cycle set4$cycle=cycle ## return return(list(raw=set2,nums=set4)) } ## run a single simulation (Figure 1) sims=par_vary(pvar=0.75, rvar=1/(3/12), tseason=c("breed"), bvar=0.05, cvar=0.5, cnvar=cf, cfvar=cf, smvar=sigma_m, swvar=sigma_w, b0var=b0, rtime='start') sim=sims$nums ## trim to final year sim=sim[which(sim$years==max(sim$years)),] ## get figure months fmonths2=c(1/12,3/12,5/12,7/12,9/12,11/12) ## determine middle of each period seasons=c(mean(c(0,winter_depart)), mean(c(winter_depart,b_start)), mean(c(b_start,b_depart)), mean(c(b_depart,winter_arrive)), mean(c(winter_arrive,1))) ## labels slabels=c('w','m','b','m','w') slabels_raw=slabels slabels=paste('italic(',slabels,')',sep='') ## plot proportion of each state ggplot(sim,aes(annual,prev),group=group2)+ ## add migration timing geom_vline(xintercept=c(winter_depart,b_start, b_depart,winter_arrive), size=0.25,linetype=2)+ ## add period annotate('text', x=seasons, y=max(sim$prev)+15, label=slabels,parse=T,size=3)+ ## add add compartment sizes geom_path(aes(colour=group2),size=0.9)+ guides(colour=F)+ ## theme theme_bw()+ theme(axis.text=element_text(size=9), axis.title=element_text(size=12), axis.text.x=element_text(angle=45,hjust=1))+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())+ theme(axis.title.x=element_text(margin=margin(t=5,r=0,b=0,l=0)))+ theme(axis.title.y=element_text(margin=margin(t=0,r=10,b=0,l=0)))+ labs(x="annual cycle", y=expression(paste("population size (",italic(S+I+L),")")))+ ## x axis scale_x_continuous( breaks=fmonths2,labels=c("1 Feb","1 Apr","1 June","1 Aug","1 Oct","1 Dec"), expand=c(0.001,0.005)) ## function to get the proportions from each list lget=function(x){ y=lapply(x,function(x) x$raw) return(y) } ## new function to get final year per parameter combination ytrim=function(x){ ## unique ydata=lapply(unique(x$upars),function(y){ ## trim set=x[which(x$upars==y),] set=set[which(set$years==max(set$years)),] return(set) }) ## trim ydata=do.call(rbind,ydata) ## return return(ydata) } ## vector of relapse evec=seq(0,1,by=0.1) ## infectious period as months rvec=1/(c(1:3)/12)[c(1,3)] ## vector of beta bvec=c(0.05,0.5) ## vector of c cvec=seq(cf,0.8,by=0.05) ## vector of winter survival swvec=c(sigma_w,sigma_b) ## vector of b0 b0vec=c(b0,b0min) ## loop for given sigma_w, transmission phenology, b0 vloop=function(slvl,trans,b){ ## outer list for varying beta oolist=list() for(k in 1:length(bvec)){ ## outer list for varying latency rate (rho) olist=list() for(j in 1:length(rvec)){ ## inner list for varying cost of infection for survival rlist=list() for(i in 1:length(cvec)){ ## run pvary, tweak other parameters (e.g., relapse at start or end) test=lapply(evec,par_vary, tseason=trans, rvar=rvec[j], bvar=bvec[k], cvar=cvec[i], cnvar=min(cvec), cfvar=min(cvec), smvar=sigma_m, swvar=swvec[slvl], b0var=b, rtime='start') test=lget(test) test=do.call(rbind,test) test=ytrim(test) rlist[[i]]=test ## print print(paste("cs = ",cvec[i],", ",i," of ",length(cvec),sep="")) } rtest=do.call(rbind,rlist) ## save outer list olist[[j]]=rtest ## print print(paste("rho = ",rvec[j],", ",j," of ",length(rvec))) } ## collapse otest=do.call(rbind,olist) ## save outer outer list oolist[[k]]=otest ## print print(paste("beta =",bvec[k],", ",k," of",length(bvec))) } ## collapse ootest=do.call(rbind,oolist) return(ootest) } ## run parameter exploration for breeding transmission, relapse at start slist=list() for(i in 1:length(swvec)){ sloop=vloop(slvl=i,trans=c("breed"),b=b0) slist[[i]]=sloop print(paste("sigma_w =",swvec[i],", ",i," of",length(swvec),"full loops finished")) } stest=do.call(rbind,slist) ## combine upars with trans stest$upars2=with(stest,paste(upars,breed_trans,mig_trans,winter_trans,timing)) ## get id values ids=stest[!duplicated(stest$upars2), c("upars2","rprop","cs","rho","sigma_m","sigma_w","beta", "breed_trans",'mig_trans','winter_trans','b0','timing')] ids=ids[order(ids$rprop),] ## get the min and max per id mm=data.frame(tapply(stest$I_prop,stest$upars2,min), tapply(stest$I_prop,stest$upars2,max), tapply(stest$N,stest$upars2,min), tapply(stest$N,stest$upars2,max), tapply(stest$years,stest$upars2,max)) names(mm)=c("p_min","p_max","N_min","N_max","ymax") mm$upars2=rownames(mm) ## merge mdata=merge(mm,ids,by="upars2",all=T) ## rho2 mdata$rho2=mdata$rho*12 mdata$rho2=(1/mdata$rho)*12 ## rprop as % mdata$rprop2=mdata$rprop*100 ## make labels mdata$rho3=paste("1/rho ==",mdata$rho2) mdata$sigma_m2=paste("sigma[m] ==",mdata$sigma_m*100,"*'%'") mdata$sigma_w2=paste("sigma[w] ==",mdata$sigma_w*100,"*'%'") mdata$beta2=paste("beta ==",mdata$beta) ## trim to one sigma value and transmission mdata2=mdata[which(mdata$breed_trans==1 & mdata$mig_trans==0 & mdata$winter_trans==0 & mdata$sigma_w==unique(mdata$sigma_w)[1]),] ## recreate Figure 3 ggplot(mdata2,aes(cs,p_max,colour=rprop2,group=rprop2))+ ## add curves geom_line()+ ## facet facet_grid(beta2~rho3,labeller= label_parsed,scales="free_y")+ ## color scale scale_colour_viridis_c(option="A",end=0.6)+ ## theme theme_bw()+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())+ labs(x=expression(paste("cost of infection for migrant survival (",italic(c[m]),")")), y="maximum steady state infection prevalence")+ guides(colour=guide_legend(title="% relapse at migration"))+ theme(legend.position="top", legend.text=element_text(size=14), legend.title=element_text(size=14))+ theme(axis.title.x=element_text(margin=margin(t=10,r=0,b=0,l=0)))+ theme(axis.title.y=element_text(margin=margin(t=0,r=10,b=0,l=0)))+ theme(axis.title=element_text(size=14), axis.text=element_text(size=11), strip.text=element_text(size=14))