################################################################################ # Course SMACH - January 2014 # a classical SEIR model. From original publication of Zadoks (1971) # without boxcar train # with weather effect on latency and infectious period # version $Date:: 2014-04-30 1#$ by $Author: brun $(v0: 2013-12-11) ################################################################################ zakoks.simpleweather.model = function (nlpd_T=40,nipd_T=200,dmfr=16,SITE0 = 1e+10,XINF0=10,weather, 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-XINF0 XLAT[sdate] <- 0 XINF[sdate] <- XINF0 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]) TMEAN= max((weather$TMIN[day]+weather$TMAX[day])/2, 0) # rapp: rate of apparition : nb of sites Latent=>Infectant rapp= min(XLAT[day] * TMEAN/nlpd_T, XLAT[day]+rocc) # rrem: rate of removal : nb of sites Infectant=>removed rrem=min(XINF[day] * TMEAN/nipd_T, 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_T=nlpd_T,nipd_T=nipd_T,dmfr=dmfr,SITE0 = SITE0))) } ################################################################################ graph_epid_s=function(out,typel="s",all=TRUE, param=TRUE,SITE0=SITE0,titre=""){ # a tip to show value=0 on the log plot sim=out$sim par(ann=TRUE,new=FALSE) 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)",main=titre) 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) duration = 50 sdate1 = 120#as.numeric(format(strptime("1997-05-25", "%Y-%m-%d", tz=""), "%j")) weather1=subset(weather_FranceWest, WEYR==1997 & idsite==39) mean( ((weather1$TMIN+weather1$TMAX)/2)[sdate1:(sdate1+duration)]) sdate2 = 120#as.numeric(format(strptime("1999-05-25", "%Y-%m-%d", tz=""), "%j")) weather2=subset(weather_FranceWest, WEYR==1999 & idsite==39) mean( ((weather2$TMIN+weather2$TMAX)/2)[sdate2:(sdate2+duration)]) sdate3 = 120#as.numeric(format(strptime("2001-05-25", "%Y-%m-%d", tz=""), "%j")) weather3=subset(weather_FranceWest, WEYR==2001 & idsite==39) mean( ((weather3$TMIN+weather3$TMAX)/2)[sdate3:(sdate3+duration)]) sdate4 = 120#as.numeric(format(strptime("2004-05-25", "%Y-%m-%d", tz=""), "%j")) weather4=subset(weather_FranceWest, WEYR==2004 & idsite==39) mean( ((weather4$TMIN+weather4$TMAX)/2)[sdate4:(sdate4+duration)]) # unlimited number of Site SITE0= 1e+10 XINF0=1/10000000 *SITE0 ################################################################################ # 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(2,2), mar=c(4,4,3,4)) out1=graph_epid_s(zakoks.simpleweather.model(nlpd_T=80,nipd_T=200,dmfr=16,SITE0 = SITE0,XINF0=XINF0, weather1, sdate = sdate1, ldate = sdate1+duration ),typel="l",SITE0= 1e+10,titre="1997") out2=graph_epid_s(zakoks.simpleweather.model(nlpd_T=80,nipd_T=200,dmfr=16,SITE0 = SITE0,XINF0=XINF0, weather2, sdate = sdate2, ldate = sdate2+duration ),typel="l",SITE0= 1e+10,titre="1999") out3=graph_epid_s(zakoks.simpleweather.model(nlpd_T=80,nipd_T=200,dmfr=16,SITE0 = SITE0,XINF0=XINF0, weather3, sdate = sdate3, ldate = sdate3+duration ),typel="l",SITE0= 1e+10,titre="2001") out4=graph_epid_s(zakoks.simpleweather.model(nlpd_T=80,nipd_T=200,dmfr=16,SITE0 = SITE0,XINF0=XINF0, weather4, sdate = sdate4, ldate = sdate4+duration ),typel="l",SITE0= 1e+10,titre="2004") ################################################################################ # compare with data fo Bancal 2008 # A) load and show data setwd("D:/data/FormationSmach/data/leaf_rust") bancal_2008_rust = read.table("leaf_rust_wheat.data.csv", sep=";", header=TRUE, dec=".", stringsAsFactors=FALSE) # select data for low Nitrogen only sdata=subset(bancal_2008_rust, Nitrogen=="low") lty_year=1:length(unique(sdata$Year)) names(lty_year)=unique(sdata$Year) table_col=c("1997"="blue", "1999"="red","2001"="green","2004"="orange") sdate = 120#as.numeric(format(strptime("1997-05-25", "%Y-%m-%d", tz=""), "%j")) dev.new() par(mfrow=c(1,1), mar=c(4,4,3,4)) plot(c(sdate,sdate+duration),c(0,100),xlab="Day of year", ylab="Disease Severity(%)", type="n") for (year in c(1997, 1999,2001,2004)){ #simulation weather=subset(weather_FranceWest, WEYR==year & idsite==39) out=zakoks.simpleweather.model(nlpd_T=80,nipd_T=200,dmfr=16,SITE0 = SITE0,XINF0=XINF0, weather, sdate = sdate, ldate = sdate+duration) lines(out$sim$day,out$sim$severity*100,col=table_col[as.character(year)],lwd=2, main=year) #observation points(sdata[sdata$Year==year,c("day","DS")],col=table_col[as.character(year)]) } legend("topleft",legend=c("sim 1997", "sim 1999","sim 2001","sim 2004","obsflagleaf 1997", "obsflagleaf 1999","obsflagleaf 2001","obsflagleaf 2004"), col=c(table_col, table_col), lty=c(1,1,1,1,NA,NA,NA,NA), pch=c(NA,NA,NA,NA,1,1,1,1), cex=0.5,ncol=2) # epidemic vs date # end of file