################################################################################ # Course SMACH - January 2014 # a classical SEIR model. From original publication of Zadoks (1971) # without boxcar train # without any weather effect # version $Date:: 2014-04-30 1#$ by $Author: brun $(v0: 2013-12-11) ################################################################################ zakoks.simple.model = function (nlpd=4,nipd=1,dmfr=16,SITE0 = 1e+10,weather=NULL, sdate = 1, ldate = 140){ XVAC <- rep(NA, ldate) # Latent XLAT <- rep(NA, ldate) # Infectious XINF <- rep(NA, ldate) XCTR <- rep(NA, ldate) # initialization of state variable XVAC[sdate] <- SITE0-1 XLAT[sdate] <- 1 XINF[sdate] <- 0 XCTR[sdate] <- 0 for (day in sdate:(ldate - 1)) { #if (day==inf_start){BOXI[day,"XINF"]=1, BOXI[sdate,"sumTEEQ"] = 1} # correction factor : feed back from total occupied sites cofr<-max(XVAC[day]/SITE0, 0) # rocc: rate of occupation : nb of sites Vacant=>Latent rocc = min(cofr * dmfr * XINF[day], XVAC[day]) # rapp: rate of apparition : nb of sites Latent=>Infectant rapp= min(XLAT[day] * 1/nlpd, XLAT[day]+rocc) # rrem: rate of removal : nb of sites Infectant=>removed rrem=min(XINF[day] * 1/nipd, XINF[day]+rapp) XVAC[day+1] <- XVAC[day]-rocc XLAT[day+1] <-XLAT[day] +rocc - rapp XINF[day+1] <- XINF[day] + rapp - rrem XCTR[day+1] <- XCTR[day] + rrem } XTO1 = XLAT+XINF+XCTR XSEV = XINF+XCTR severity=XSEV/(XLAT+XINF+XCTR+XVAC) return(list(sim=data.frame(day = sdate:ldate, XVAC = XVAC[sdate:ldate], XLAT = XLAT[sdate:ldate], XINF = XINF[sdate:ldate],XCTR = XCTR[sdate:ldate],XTO1=XTO1[sdate:ldate], XSEV= XSEV[sdate:ldate], severity=severity[sdate:ldate]), param=c(nlpd=nlpd,nipd=nipd,dmfr=dmfr,SITE0 = SITE0))) } ################################################################################ graph_epid_s=function(out,typel="s",all=TRUE, param=TRUE,SITE0=SITE0){ # a tip to show value=0 on the log plot sim=out$sim plot(sim$day,sim$XSEV, type=typel,lty=1, lwd=3, xlab="Time after first infection (day)", ylim=c(0,SITE0),ylab="X (nb of infectious sites)") if(all){ lines(sim$day,sim$XLAT, type=typel,lty=1, lwd=2, col="orange") lines(sim$day,sim$XINF, type=typel,lty=1,lwd=2, col="red") lines(sim$day,sim$XCTR, type=typel,lty=1,lwd=2, col="grey") } if(param){ text(seq(30,100,length.out=length(out$param)),max(sim$XSEV),names(out$param),cex=0.75) text(seq(30,100,length.out=length(out$param)),0.95*max(sim$XSEV),(out$param),cex=0.75) } par(ann=FALSE,new=TRUE) plot.default(sim$day,sim$severity,axes=FALSE,type="l",ylim=c(0,1),ylab='',lty=3, lwd=3,col="darkgrey") axis(side=4) mtext("severity (%)", side=4,line=2) mtext(paste("AUDPC = ",round(sum(sim$severity,1))),side=3,line=0,cex=0.8) legend("topleft",legend=c("XSEV","XLAT","XINF","XCTR","severity"),lty=c(1,1,2,1,3), lwd=c(2,2,2,2,3),col=c("black","orange","red","grey","darkgrey"), cex=0.6, ncol = 2) return(invisible(out)) } ################################################################################ library(ZeBook) # unlimited number of Site P.recondita = brown rust =1e+10 sites for LAI=4 with 320stems/m2 SITE0= 1e+10 sdate = 1 ldate = 150 out=zakoks.simple.model(nlpd=16,nipd=30,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = ldate ) ################################################################################ # Latent period varies from 8-14 days in the field with temperatures of 10-25 C. # With temperatures of 5-10 C latent periods of several weeks are common. # typical infectious period for leaf rust about 30 days par(mfrow=c(3,1), mar=c(4,4,2,4)) graph_epid_s(zakoks.simple.model(nlpd=8,nipd=30,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = ldate ),typel="l",SITE0=SITE0) graph_epid_s(zakoks.simple.model(nlpd=16,nipd=30,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = ldate ),typel="l",SITE0=SITE0) graph_epid_s(zakoks.simple.model(nlpd=24,nipd=30,dmfr=16,SITE0 = SITE0, weather, sdate = sdate, ldate = ldate ),typel="l",SITE0=SITE0) ################################################################################ # end of file