Initial commit

.gitignore

+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+*.Rproj

Scripts/dynamic.R

+#header #################################################################################
+#'Dynamic_model.R'
+
+#Title: MADRA
+#Project 1600001916: pid
+#Client: client
+#Author: <Eduardo> <Costa>, Wageningen Bioveterinary Research
+
+#Description: Multi-directional antimicrobial resistance modeling
+
+#Start date: 9/1/2020
+#Last Update: {6:date}
+
+#R version: r.version
+#Scriptversion: version
+
+#Dependencies
+#<-Downstream
+#->Upstream
+
+#Input:
+#-
+
+#Output:
+#-
+
+#Peer reviewer(s)
+
+#Please ensure directories are relative. Please comment code sufficiently.
+
+#Script start#############################################################################
+
+model<-function(iteraction,times,flocks){
+
+
+simtable=8
+
+#Organizing the output directory
+
+wd <- getwd()
+
+data.labels <- paste("data",1:iteraction, ".txt", sep="")
+
+
+layout.matrix <- matrix(c(1, 2, 3, 4,0,0), byrow=T,nrow = 3, ncol = 2)
+
+layout(mat = layout.matrix,
+       heights = c(2, 2,0.5), # Heights of the 3 rows
+       widths = c(2, 2)) # Widths of the two columns
+
+
+
+for(j in 1:simtable){ # looping Scenarios
+ # All inputs are stored in this file
+
+  source("inputs.R")
+   
+  wd.dados <- paste(wd, "/output_scenario",j,sep="")
+  dir.create( paste(wd, "/output_scenario",j, sep="") )
+  
+  par(mar=c(5,5,4,3))
+  plot(1:times,seq(0,1,length.out = (times)),type="n", 
+       ylab="Prevalence", xlab="Time (Weeks)",xlim=c(1,times),cex=1.2, cex.main=2, cex.lab=2, cex.axis=2)
+ 
+   #main=paste(j,":","f_hu=",fade_h[j],"&","f_ch=",fade_ch[j],"&","b_ch=",beta_ch[j]),
+  
+  
+  for(iter in 1:iteraction){ #Stochasticity at time t
+    tempo_inicial1 = proc.time()
+    cat("\n\n iteraction: ", iter, "\n")
+    
+    source("inputs.R")
+    rdu<-function(n,k) sample(1:k,n,replace=T)
+  
+    H[1]<-prev1
+    Fa[1]<-prev2
+    time_flock[[1]]<-rdu(flocks,6)
+    
+    
+    
+      r1[1]<-  (H[1]*(1-H[1]))*col_proc
+      r2[1]<-  (Fa[1])*(1-H[1])*expo_fa[j]
+      r3[1]<-(1 - (1 + ((dose_chi[1] + dose_veg[1])/alfa))^-gama)*non_veg*prob_consu[j]*prev_chi[1]
+      r4[1]<-(1 - (1 + ((dose_veg[1])/alfa))^-gama)*non_veg*prob_consu[j]*prev_veg[1]
+      r5[1]<- (1-exp(-1/fade_h[j])) # decolonization probability in one week  
+      
+      p1[1]<-(r1[1]+r2[1]+r3[1]+r4[1])-(r1[1]*r2[1])-(r1[1]*r3[1])-(r1[1]*r4[1])-(r2[1]*r3[1])-(r2[1]*r4[1])-(r3[1]*r4[1])
+      +(r1[1]*r2[1]*r3[1])+(r1[1]*r2[1]*r4[1])+(r1[1]*r3[1]*r4[1])+(r2[1]*r3[1]*r4[1])-(r1[1]*r2[1]*r3[1]*r4[1]) 
+      p2[1]<-r5[1]
+
+      r6[1]<-(1-Fa[1])*(H[1])*col_proc2
+      r7[1]<-(1-Fa[1])*(C[1])*rate1[j]
+      r8[1]<-r4[1]
+      
+      p3[1]<-(r6[1]+r7[1]+r3[1]+r4[1])-(r6[1]*r7[1])-(r6[1]*r3[1])-(r6[1]*r4[1])-(r7[1]*r3[1])-(r7[1]*r4[1])-(r3[1]*r4[1])
+      +(r6[1]*r7[1]*r3[1])+(r6[1]*r7[1]*r4[1])+(r6[1]*r3[1]*r4[1])+(r7[1]*r3[1]*r4[1])-(r6[1]*r7[1]*r3[1]*r4[1]) 
+      p4[1]<- r5[1]    
+
+      
+      r9[[1]]<-  (1-exp(-beta_ch[j]*time_flock[[1]]))*prev3
+      r10[[1]]<-  (Fa[1])*(1-prev3)*rate2[j]
+      r11[[1]]<- (1-exp(-(1/fade_ch[j])))
+      p5[[1]]<-  (r9[[1]]+r10[[1]])-(r9[[1]]*r10[[1]])
+      p6[[1]]<-  r11[[1]]
+     
+      flock[[1]]<- (1-prev3)*p5[[1]]+prev3*(1-p6[[1]]) 
+      C[1]<-mean(flock[[1]])
+      media[1]<-mean(flock[[1]][time_flock[[1]]>5])
+
+
+    Nenv1[1]<-(  media[1]*aext1*10^(conc1*afec1+log(10)*sig_conc1/2) + 
+                   afec1*pfec1*wfec1*(media[1])*10^(conc1+log(10)*sig_conc1/2) )/((benv1+cenv1)-(benv1*cenv1))
+    
+    Next1[1]<-(1-aext1)*(1-cext1)*media[1]*10^(conc1*afec1+log(10)*sig_conc1/2 )+
+      (benv1*(Nenv1[1]))+
+      ((1-afec1)*pfec1*wfec1*(media[1])*10^(conc1+log(10)*sig_conc1/2))
+    
+    #2 step in slaughterhouse: DEFEATHERING
+    
+    Nenv2[1]<-(  aext2*(Next1[1]) + afec2*wfec2*pfec2*(media[1])*10^(conc2+log(10)*sig_conc2/2) )/((benv2+cenv2)-(benv2*cenv2))
+    
+    Next2[1]<-(1-aext2)*(1-cext2)*(Next1[1])+
+      (benv2*(Nenv2[1]))+
+      (1-afec2)*wfec2*pfec2*(media[1])*10^(conc2+log(10)*sig_conc2/2)
+    
+    
+    #3 step in slaughterhouse: EVISCERATION
+    
+    
+    Nenv3[1]<-(  aext3*(Next2[1]) + afec3*pfec3*(media[1])*wfec3*10^(conc3+log(10)*sig_conc3/2) )/((benv3+cenv3)-(benv3*cenv3))
+    
+    Next3[1]<-(1-aext3)*(1-cext3)*(Next2[1])+
+      (benv3*(Nenv3[1]))+
+      (1-afec3)*wfec3*pfec3*(media[1])*10^(conc3+log(10)*sig_conc3/2)
+    
+    
+    #5 step in slaughterhouse: CHILLING
+    
+    
+    Nenv5[1]<- ( aext5*(Next3[1]) )/((benv5+cenv5)-(benv5*cenv5))
+    
+    Next5[1]<-(1-aext5)*(1-cext5)*(Next3[1])+
+      (benv5*(Nenv5[1]))
+    
+    
+    if (Next5[1]==0){
+      
+      dose_chi[1]<-0
+      dose_veg[1]<-0
+      
+    }else{
+      
+      
+      a[1]<-log10(Next5[1]/car_me) #log10(cfu/g)
+      prev_slaug[1]<-(media[1]+(1-exp(-Next5[1]/car_me)))-(media[1]*(1-exp(-Next5[1]/car_me)))
+      
+      
+      #Sub-module Pre-retail reduction#  Evers`s model#
+      #################################################
+      
+      
+      raw[1]<-(a[1]) #log10(cfu/g) concentration at the time t after retail phase. Notice no reduction is happenig
+      
+      
+      
+      #Sub-module growth at consumer# Evers`s model#
+      ########################################################
+      
+      tgen<-(tgen_min/((Tst-Tmin)/(Topt-Tmin))^2)
+      g<-((log(2)/tgen)*tst)/log(10)   #Multiplication factor at 18 C and for 10.8 minutes
+      
+      
+      raw_af[1]<-a[1]+g  # log10(cfu/g) concentration in time t after storage at 18 C and for 10.8 minutes
+      
+      
+      
+      #Sub-module cross-contamination at consumer# Evers`s model#
+      ###########################################################
+      
+      raw[1]<- (10^raw_af[1])*consu_ch*(1-frac_cross)
+      
+      raw_cross[1]<-(10^raw_af[1])*consu_ch * (1-cross) *(frac_cross)*prev_slaug[1]  #CFU on 1 portion
+      
+      
+      veg_cross[1]<-(10^raw_af[1])*consu_ch * (cross)*frac_cross * prev_slaug[1]      #CFU on 1 portion
+      
+      
+      
+      
+      
+      #Sub-module inactivation at consumer# Evers`s model#
+      ########################################################
+      
+      R<-the*10^((Tend-Tref)/z)/Dref
+      
+      #Reduction factor at 71.5 C and 15 minutes
+      
+      
+      ready[1]<-log10(raw_cross[1]/consu_ch)+ R    #Final concentration log10(cfu/g) in 1 serving
+      
+      ready2[1]<-(log10(raw[1]/consu_ch)+ R)  #Final concentration log10(cfu/g) in 1 serving without salad
+      
+      
+      
+      
+      
+      #Sub-module exposure cfu#
+      ###########################
+      dose_chi[1]<-( (10^ready[1])+(10^ready2[1]) )*consu_ch     #total dose/serving (cfu)
+      
+      dose_veg[1]<-(veg_cross[1])*1.8  #total dose/serving (cfu)
+      
+      prev_chi[1]<-1-exp(-10^dose_chi[1])
+      
+      prev_veg[1]<-1-exp(-veg_cross[1]*1.8)
+      
+    } #end of if steatment
+    
+    
+    
+    
+    
+    for (t in 2:times){ #Time running in times 
+      tempo_inicial2 = proc.time()
+      cat("\n\n day: ", t, "\n")
+      
+      
+      ####################################
+      #Module Human and animal population#
+      ####################################
+      
+      
+      
+      #Sub-module-Open community population#
+      #####################################
+      
+      
+      #Frac of open community infec by open community + Frac of infec by farmers +Frac. of open community infec by food. 
+      
+      r1[t]<-  (H[t-1]*(1-H[t-1]))*col_proc
+      r2[t]<-  (Fa[t-1])*(1-H[t-1])*expo_fa[j]
+      r3[t]<-(1 - (1 + ((dose_chi[t-1] + dose_veg[t-1])/alfa))^-gama)*non_veg*prob_consu[j]*prev_chi[t-1]
+      r4[t]<-(1 - (1 + ((dose_veg[t-1])/alfa))^-gama)*non_veg*prob_consu[j]*prev_veg[t-1]
+      r5[t]<- (1-exp(-1/fade_h[j])) # decolonization probability in one week  
+      
+      p1[t]<-(r1[t]+r2[t]+r3[t]+r4[t])-(r1[t]*r2[t])-(r1[t]*r3[t])-(r1[t]*r4[t])-(r2[t]*r3[t])-(r2[t]*r4[t])-(r3[t]*r4[t])
+      +(r1[t]*r2[t]*r3[t])+(r1[t]*r2[t]*r4[t])+(r1[t]*r3[t]*r4[t])+(r2[t]*r3[t]*r4[t])-(r1[t]*r2[t]*r3[t]*r4[t]) 
+      p2[t]<-r5[t]
+      
+      
+      H[t]<-  (1-H[t-1])*p1[t]+ H[t-1]*(1-p2[t])
+      
+      
+      #Sub-module-farmers population#
+      ###############################
+       #Frac of farmers infec by open comunity+ Frac of farmers infec by chickens + Frac of farmers infec by food
+
+	      r6[t]<-(1-Fa[t-1])*(H[t-1])*col_proc2
+        
+        r7[t]<-(1-Fa[t-1])*(C[t-1])*rate1[j]
+        
+        r8[t]<-r5[t]
+
+        p3[t]<-(r6[t]+r7[t]+r3[t]+r4[t])-(r6[t]*r7[t])-(r6[t]*r3[t])-(r6[t]*r4[t])-(r7[t]*r3[t])-(r7[t]*r4[t])-(r3[t]*r4[t])
+        +(r6[t]*r7[t]*r3[t])+(r6[t]*r7[t]*r4[t])+(r6[t]*r3[t]*r4[t])+(r7[t]*r3[t]*r4[t])-(r6[t]*r7[t]*r3[t]*r4[t]) 
+
+	      p4[t]<- r8[t]    
+	
+      	Fa[t]<-  (1-Fa[t-1])*p3[t]+ Fa[t-1]*(1-p4[t])
+      
+      #Sub-module-chicken Population#
+      ###############################
+      
+      
+      #Frac. of chikens infect by chickens + #Frac. of chickens infect by farmers - Carries status self clearance (times)
+      
+      	time_flock[[t]]<-ifelse (time_flock[[t-1]]<6,time_flock[[t-1]]+1,1)
+      
+      	r9[[t]]<-ifelse(time_flock[[t]]==1,(1-exp(-beta_ch[j]*prev3)),(1-exp(-beta_ch[j]*flock[[t-1]])) )
+      	r10[[t]]<-ifelse(time_flock[[t]]!=1,(Fa[t-1])*(1-flock[[t-1]])*rate2[j],  (Fa[t-1])*(1-prev3)*rate2[j])
+      	r11[[t]]<-ifelse(time_flock[[t]]!=1,(1-exp(-(1/fade_ch[j]))), (1-exp(-(1/fade_ch[j]))))
+      	p5[[t]]<-ifelse(time_flock[[t]]!=1,(r9[[t]]+r10[[t]])-(r9[[t]]*r10[[t]]), (r9[[t]]+r10[[t]])-(r9[[t]]*r10[[t]]))
+      	p6[[t]]<-ifelse(time_flock[[t]]!=1,r11[[t]], r11[[t]])
+      	flock[[t]]<-ifelse(time_flock[[t]]!=1,(1-flock[[t-1]])*p5[[t]]+flock[[t-1]]*(1-p6[[t]]), (1-prev3)*p5[[t]]+prev3*(1-p6[[t]])  )
+  
+      	
+      C[t]<-mean(flock[[t]])
+      media[t]<-mean(flock[[t]][time_flock[[t]]>5])     
+
+      
+      ###########################
+      #Module food contamination#
+      ###########################
+      
+      #Sub-module slaughterhouse# Nauta`s model#
+      #########################################
+      
+      ##Expected value
+      #1 step in slaughterhouse: SCALDING
+      
+      Nenv1[t]<-(  media[t]*aext1*10^(conc1*afec1+log(10)*sig_conc1/2) + 
+                     afec1*pfec1*wfec1*(media[t])*10^(conc1+log(10)*sig_conc1/2) ) /((benv1+cenv1)-(benv1*cenv1))
+      
+      Next1[t]<-(1-aext1)*(1-cext1)*media[t]*10^(conc1*afec1+log(10)*sig_conc1/2 )+
+        (benv1*(Nenv1[t]))+
+        ((1-afec1)*pfec1*wfec1*(media[t])*10^(conc1+log(10)*sig_conc1/2))
+      
+      #2 step in slaughterhouse: DEFEATHERING
+      
+      Nenv2[t]<-(  aext2*(Next1[t]) + afec2*wfec2*pfec2*(media[t])*10^(conc2+log(10)*sig_conc2/2) )/((benv2+cenv2)-(benv2*cenv2))
+      
+      Next2[t]<-(1-aext2)*(1-cext2)*(Next1[t])+
+        (benv2*(Nenv2[t]))+
+        (1-afec2)*wfec2*pfec2*(media[t])*10^(conc2+log(10)*sig_conc2/2)
+      
+      
+      #3 step in slaughterhouse: EVISCERATION
+      
+      
+      Nenv3[t]<-(  aext3*(Next2[t]) + afec3*pfec3*(media[t])*wfec3*10^(conc3+log(10)*sig_conc3/2) )/((benv3+cenv3)-(benv3*cenv3))
+      
+      Next3[t]<-(1-aext3)*(1-cext3)*(Next2[t])+
+        (benv3*(Nenv3[t]))+
+        (1-afec3)*wfec3*pfec3*(media[t])*10^(conc3+log(10)*sig_conc3/2)
+      
+     
+      #4 step in slaughterhouse: CHILLING
+      
+      
+      Nenv5[t]<- ( aext5*(Next3[t]) )/((benv5+cenv5)-(benv5*cenv5))
+      
+      Next5[t]<-(1-aext5)*(1-cext5)*(Next3[t])+
+        (benv5*(Nenv5[t]))
+      
+      
+      if (Next5[t]==0){
+        
+        dose_chi[t]<-0
+        dose_veg[t]<-0
+        
+      }else{
+        
+        
+        a[t]<-log10(Next5[t]/car_me) #log10(cfu/g)
+        
+        prev_slaug[t]<-((1-exp(-Next5[t]/car_me))+media[t])-((1-exp(-Next5[t]/car_me))*media[t])
+        
+        #Sub-module Pre-retail reduction#  Evers`s model#
+        #################################################
+        
+        
+        #(a[t]) #log10(cfu/g) concentration at the time t after retail phase. Notice no reduction is happenig
+        
+        
+        
+        #Sub-module growth at consumer# Evers`s model#
+        ########################################################
+        
+        tgen<-(tgen_min/((Tst-Tmin)/(Topt-Tmin))^2)
+        g<-((log(2)/tgen)*tst)/log(10)   #Multiplication factor at 18 C and for 10.8 minutes
+        
+        
+        raw_af[t]<-a[t]+g  # log10(cfu/g) concentration in time t after storage at 18 C and for 10.8 minutes
+        
+        
+        
+        #Sub-module cross-contamination at consumer# Evers`s model#
+        ###########################################################
+        
+        raw[t]<- (10^raw_af[t])*consu_ch*(1-frac_cross)
+        
+        raw_cross[t]<-(10^raw_af[t])*consu_ch * (1-cross) *(frac_cross)*prev_slaug[t]  #CFU on 1 portion
+        
+        
+        veg_cross[t]<-(10^raw_af[t])*consu_ch * (cross)*frac_cross * prev_slaug[t]      #CFU on 1 portion
+        
+        
+        
+        
+        #Sub-module inactivation at consumer# Evers`s model#
+        ########################################################
+        
+        R<-the*10^((Tend-Tref)/z)/Dref
+        
+        #Reduction factor at 71.5 C and 15 minutes
+        
+        
+        ready[t]<-log10(raw_cross[t]/consu_ch)+ R    #Final concentration log10(cfu/g) in 1 serving with salad
+        
+        ready2[t]<-(log10(raw[t]/consu_ch)+ R)  #Final concentration log10(cfu/g) in 1 serving without salad
+        
+        
+        #Sub-module exposure cfu#
+        #########################
+        
+        dose_chi[t]<-( (10^ready[t])+(10^ready2[t]) )*consu_ch     #total dose/serving (cfu)
+        
+        dose_veg[t]<-(veg_cross[t])*1.8  #total dose/serving (cfu)
+        
+        prev_chi[t]<-1-exp(-10^dose_chi[t])
+        
+        prev_veg[t]<-1-exp(-veg_cross[t]*1.8)
+        
+      } #end of if steatment
+      
+      
+      #Setting up outputs#
+      ####################
+      
+      
+    } # End of the looping (time t)
+    
+    
+    data<- cbind(H,media,C,Fa)
+    
+    
+    data.local  <- file.path(wd.dados, data.labels[iter])
+    
+    
+    write.table(data, file = data.local, row.names=FALSE, col.names=TRUE, append=FALSE)
+    
+    lines(1:times, H, col = "red")
+    lines(1:times, C, col = "green")
+    lines(1:times, Fa, col = "blue")
+    lines(1:times, media, col = "yellow")
+    
+    media_H[iter]<-c(H[times])
+    
+  } # End of the looping (stochasticity on the time t)
+  
+  media_H2[j]<-mean(media_H)  
+  
+  
+}#End of scenarios
+
+par(xpd=NA)
+legend(-65,-0.5, legend=c("Ope. Commu.", "Chicken","Farmer","Slaug_prev"),
+       horiz=TRUE,col=c("red","green", "blue","yellow"), lty=1:1, cex=1.2)
+
+
+saveRDS(media_H2,"cena.rds")
+
+
+
+cena<-readRDS("cena.rds")
+
+tornado<-function(){
+  par(mfrow=c(1,1))
+  base<-cena[1]
+nome<-c("fade_h=14", "fade_h=8", "fade_ch=0.5", "beta_ch=5", "rate=0.1","Very_low","No_ch")
+
+dado<-cbind.data.frame(cena[2:8],nome)
+
+sorted<-dado[order(abs(log(cena[2:8]/base))),]
+
+par(mar=c(5,8,1,2))
+barplot(log(sorted$cena/base),  horiz = T, xlim=c(-2.5,0.1),
+        names.arg = sorted$nome,col = "darkgray",las=1,ylab="Scenarios",xlab="Log reduction",
+        ann = FALSE)
+
+mtext(side = 1, text = "log reduction", line = 3)
+
+mtext(side = 2, text = "Scenarios", line = 6.5)
+}
+tornado()
+
+}#end function
+

main.R

+#header
+
+library(here)
+
+here("output")
+
+dir.create("Output",showWarnings = F)
+dir.create("Figures",showWarnings = F)
+
+source(file.path("Scripts","dynamic.R"))
+
+rm(list=ls())