diff --git a/GeneralizedImmunityExp.R b/GeneralizedImmunityExp.R
deleted file mode 100644
index 1392ec1..0000000
--- a/GeneralizedImmunityExp.R
+++ /dev/null
@@ -1,1238 +0,0 @@
-require(deSolve)
-library(tidyverse)
-args <- commandArgs(trailingOnly = TRUE)
-print(args)
-paramFile<-args[1]
-No<-as.integer(args[2])
-resultsFolder <-args[3]
-jobid<-as.integer(args[4])
-#paramFile<-"paramList.csv"
-#No<-1
-#resultsFolder<-"./"
-outFile<-paste(resultsFolder,"/",strsplit(paramFile,split="[.]")[[1]][1],"_results.csv",sep="")
-
-SEAI_p<-function(t,y,p){
-  
-  with(as.list(c(y,p)),{
-    # total host population size
-    #print(t)
-    N = NU+ND+SU+SD+AU+AD
-    # total parasite population size
-    #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD
-    #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD
-    #P = PW+PR
-    
-    # parasite/host ratio
-    rNU <-  (PW_NU+PR_NU)/(NU+1e-6) 
-    rND <-  (PW_ND+PR_ND)/(ND+1e-6) 
-    rSU <-  (PW_SU+PR_SU)/(SU+1e-6) 
-    rSD <-  (PW_SD+PR_SD)/(SD+1e-6) 
-    rAU <-  (PW_AU+PR_AU)/(AU+1e-6) 
-    rAD <-  (PW_AD+PW_AD)/(AD+1e-6) 
-    
-    # prevalence
-    I_NU = 1-exp(-rNU)
-    I_ND = 1-exp(-rND)
-    I_SU = 1-exp(-rSU)
-    I_SD = 1-exp(-rSD)
-    I_AU = 1-exp(-rAU)
-    I_AD = 1-exp(-rAD)
-    
-    # proportion of hosts that have no parasites from W or 
-    p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) 
-    p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) 
-    p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) 
-    p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) 
-    p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) 
-    p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6))) 
-    
-    # contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite
-    ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2]))
-    ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2]))
-    ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2]))
-    ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2]))
-    ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2]))
-    ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2]))
-    
-    # converted contribution of W or R with mixed infection cos----
-    h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
-    h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost))
-    h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
-    h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
-    h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
-    h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost))
-    
-    #transmissibility is reduced to the number of infected hosts, 
-    #receiving hosts needs to have open "resource room" for accepting new strains
-    B_PW = b*(
-      h_NU[1] +
-        h_ND[1] + 
-        h_SU[1] +
-        h_SD[1] +
-        h_AU[1] +
-        h_AD[1] 
-    )
-    
-    
-    B_PR = b*(
-      h_NU[2]+
-        h_ND[2]+ 
-        h_SU[2]+
-        h_SD[2]+
-        h_AU[2]+
-        h_AD[2]
-    )
-    
-    
-    
-    # per capita biting rate for hosts who are not treated
-    Bu = (B_PW+B_PR)
-    # per capita biting rate for hosts who are drug treated
-    Bd = B_PR
-    #print(B)
-    
-    # pop size +1 is to prevent N = 0 for undefined scenario
-    fracN = P_acc_N/(NU+ND+1e-6)
-    fracS = P_acc_S/(SU+SD+1e-6)
-    fracA = P_acc_A/(AU+AD+1e-6)
-    
-    # probability that the infected strain has not been seen by the adaptive immune system
-    #infN = (1-1/nstrains)^fracN
-    infN = 1 #for generalized immunity model, infectivity values are fixed over time
-    #infS = (1-1/nstrains)^fracS
-    infS = 1#red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5
-    #infA = (1-1/nstrains)^fracA
-    infA = 1
-    #print(infN)
-    #naive class getting infected, and treated
-    NU2ND = (Bu*(1-rNU/capacity)*infN*d1)
-    
-    #naive U class getting enough # of infections to move to SU
-    NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN))))
-    #naive U class die from infections
-    NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU_gen_N*muNU)
-    
-    #naive D class getting enough # of infections to move to SD
-    ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN))))
-    
-    #naive D class die from infections
-    NDInfDie = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND)
-    
-    #SU class getting infected, show symptoms, and treated
-    SU2SD = (Bu*(1-rSU/capacity)*infS*d2)
-    
-    #SU class getting enough # of infections to move to AU
-    SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS))))
-    
-    #SD class getting enough # of infections to move to AD
-    SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS))))
-    
-    #AU class getting infected, show symptoms, and treated
-    AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)
-    
-    muWU_gen_N = muWU*exp(-red*P_acc_N)+1
-    muWD_gen_N = muWD*exp(-red*P_acc_N)+1
-    muRU_gen_N = muRU*exp(-red*P_acc_N)+1
-    muRD_gen_N = muRD*exp(-red*P_acc_N)+1
-    
-    muWU_gen_S = muWU*exp(-red*P_acc_S)+1
-    muWD_gen_S = muWD*exp(-red*P_acc_S)+1
-    muRU_gen_S = muRU*exp(-red*P_acc_S)+1
-    muRD_gen_S = muRD*exp(-red*P_acc_S)+1
-    
-    muWU_gen_A = muWU*exp(-red*P_acc_A)+1
-    muWD_gen_A = muWD*exp(-red*P_acc_A)+1
-    muRU_gen_A = muRU*exp(-red*P_acc_A)+1
-    muRD_gen_A = muRD*exp(-red*P_acc_A)+1
-    
-    # Parasite pop ODE------
-    ## sensitive-------
-    dPW_NU.dt = B_PW*NU*infN*(1-rNU/capacity) + #new infections
-      (PW_ND/Tau)- # drug loose effectiveness
-      (PW_NU*NU2ND) -  #naive class getting infected, and treated
-      (PW_NU*NU2SU) - #naive class getting enough # of infections to move to S
-      (PW_NU*NUInfDie)- #naive class die from infections
-      (PW_NU*att)- #normal death rate
-      (PW_NU*muWU_gen_N)#clearance of infection
-    
-    
-    #print(dNU.dt)
-    
-    dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) + #new infections
-      (PW_NU*NU2ND) - #naive class getting infected, and treated
-      (PW_ND/Tau) - # drug loose effectiveness
-      (PW_ND*ND2SD) -#naive class getting enough # of infections to move to S
-      (PW_ND*NDInfDie) -#naive class die from infections
-      (PW_ND*att)- #normal death rate
-      (PW_ND*muWD_gen_N)#clearance from drug
-    
-    #print(dND.dt)
-    
-    dPW_SU.dt = B_PW*SU*infS*(1-rSU/capacity) + #new infections
-      (PW_NU*NU2SU) +#naive class getting enough # of infections to move to S
-      (PW_SD/Tau)-  # drug loose effectiveness
-      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PW_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
-      (PW_SU*att)- #normal death rate
-      (PW_SU*muWU_gen_S)#clearance of infection
-    
-    dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections
-      (PW_ND*ND2SD)+ #naive class getting enough # of infections to move to S
-      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PW_SD/Tau) - # drug loose effectiveness
-      (PW_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
-      (PW_SD*att)- #normal death rate
-      (PW_SD*muWD_gen_S)#clearance from drug
-    
-    dPW_AU.dt = B_PW*AU*infA*(1-rAU/capacity) + #new infections
-      (PW_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
-      (PW_AD/Tau)-  # drug loose effectiveness
-      (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PW_AU*att)- #normal death rate
-      (PW_AU*muWU_gen_A)#clearance from infection
-    
-    dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) +  #new infections
-      (PW_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
-      (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PW_AD/Tau)-  # drug loose effectiveness
-      (PW_AD*att)- #normal death rate
-      (PW_AD*muWD_gen_A)#clearance from drug
-    
-    ## resistant-------
-    dPR_NU.dt = B_PR*NU*infN*(1-rNU/capacity) + #new infections
-      (PR_ND/Tau)- # drug loose effectiveness
-      (PR_NU*NU2ND) -  #naive class getting infected, and treated
-      (PR_NU*NU2SU) - #naive class getting enough # of infections to move to S
-      (PR_NU*NUInfDie)- #naive class die from infections
-      (PR_NU*att)- #normal death rate
-      (PR_NU*muRU_gen_N)#clearance of infection
-    
-    #print(dNU.dt)
-    
-    dPR_ND.dt = B_PR*ND*infN*(1-rND/capacity) + #new infections
-      (PR_NU*NU2ND) - #naive class getting infected, and treated
-      (PR_ND/Tau) - # drug loose effectiveness
-      (PR_ND*ND2SD) -#naive class getting enough # of infections to move to S
-      (PR_ND*NDInfDie) -#naive class die from infections
-      (PR_ND*att)- #normal death rate
-      (PR_ND*muRD_gen_N)#clearance from drug
-    
-    #print(dND.dt)
-    
-    dPR_SU.dt = B_PR*SU*infS*(1-rSU/capacity) + #new infections
-      (PR_NU*NU2SU) +#naive class getting enough # of infections to move to S
-      (PR_SD/Tau)-  # drug loose effectiveness
-      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PR_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
-      (PR_SU*att)- #normal death rate
-      (PR_SU*muRU_gen_S)#clearance of infection
-    
-    dPR_SD.dt = B_PR*SD*infS*(1-rSD/capacity) + #new infections
-      (PR_ND*ND2SD)+ #naive class getting enough # of infections to move to S
-      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PR_SD/Tau) - # drug loose effectiveness
-      (PR_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
-      (PR_SD*att)- #normal death rate
-      (PR_SD*muRD_gen_S)#clearance from drug
-    
-    dPR_AU.dt = B_PR*AU*infA*(1-rAU/capacity) + #new infections
-      (PR_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
-      (PR_AD/Tau)-  # drug loose effectiveness
-      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PR_AU*att)- #normal death rate
-      (PR_AU*muRU_gen_A)#clearance from infection
-    
-    dPR_AD.dt = B_PR*AD*infA*(1-rAD/capacity) + #new infections
-      (PR_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
-      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PR_AD/Tau)-  # drug loose effectiveness
-      (PR_AD*att)- #normal death rate
-      (PR_AD*muRD_gen_A)#clearance from drug
-    
-    # Generalized immunity ODE------
-    dNU.dt = birth + 
-      (ND/Tau)- # drug loose effectiveness
-      NU*NU2ND -  #naive class getting infected, and treated
-      NU*NU2SU - #naive class getting enough # of infections to move to S
-      NU*NUInfDie- #naive class die from infections
-      (NU*att)#normal death rate
-    
-    
-    #print(dNU.dt)
-    
-    dND.dt = NU*NU2ND - #naive class getting infected, and treated
-      (ND/Tau) - # drug loose effectiveness
-      ND*ND2SD -#naive class getting enough # of infections to move to S
-      ND*NDInfDie -#naive class die from infections
-      (ND*att) #normal death rate
-    #print(dND.dt)
-    
-    dSU.dt =  NU*NU2SU +#naive class getting enough # of infections to move to S
-      (SD/Tau)-  # drug loose effectiveness
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      SU*SU2AU -#symptomatic class getting enough # of infections to move to AU
-      (SU*att) #normal death rate
-    
-    dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      (SD/Tau) - # drug loose effectiveness
-      SD*SD2AD- #symptomatic class getting enough # of infections to move to AD
-      (SD*att) #normal death rate
-    
-    dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS
-      (AD/Tau)-  # drug loose effectiveness
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AU*att) #normal death rate
-    
-    dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AD/Tau)-  # drug loose effectiveness
-      (AD*att) #normal death rate
-    
-    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
-    
-    
-    # track the mean number of times infected per immunity class------
-    dP_acc_N.dt = NU*Bu*(1-rNU/capacity)*(1-infN) + ND*Bd*(1-rND/capacity)*(1-infN)+
-      NU*Bu*(1-rNU/capacity)*infN*muWU_gen_N + ND*Bd*(1-rND/capacity)*infN*muRD_gen_N- # total number of infections received
-      (NU*(att+NUInfDie))*fracN- #when naive U class die how many infections deducted
-      (ND*(att+NDInfDie))*fracN- #when naive D class die many infections deducted
-      (NU*NU2SU+ND*ND2SD)*fracN #when naive class leave to S
-    
-    dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S
-      SU*Bu*(1-rSU/capacity)*(infS*muWU_gen_S+(1-infS)) + 
-      SD*Bd*(1-rSD/capacity)*(infS*muRD_gen_S+(1-infS))- #total number of new infections
-      att*P_acc_S- # normal death
-      (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic
-    
-    dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic
-      AU*Bu*(1-rAU/capacity)*(infA*muWU_gen_A+(1-infA)) + 
-      AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections
-      (ImmLost+att)*P_acc_A # normal death
-    
-    a <- 365/(2*pi)
-    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
-    #db.dt =0
-    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
-    #out1<- list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
-    #              dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
-    #              dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
-    #              dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt))
-    #print(out1)
-    
-    return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
-                  dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
-                  dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
-                  dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt)))
-  })
-}
-
-SEAI_resOnly<-function(t,y,p){
-  
-  with(as.list(c(y,p)),{
-    # total host population size
-    N = NU+ND+SU+SD+AU+AD
-    # total parasite population size
-    #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD
-    #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD
-    #P = PW+PR
-    
-    # parasite/host ratio
-    rNU <-  (PW_NU+PR_NU)/(NU+1e-6) 
-    rND <-  (PW_ND+PR_ND)/(ND+1e-6) 
-    rSU <-  (PW_SU+PR_SU)/(SU+1e-6) 
-    rSD <-  (PW_SD+PR_SD)/(SD+1e-6) 
-    rAU <-  (PW_AU+PR_AU)/(AU+1e-6) 
-    rAD <-  (PW_AD+PW_AD)/(AD+1e-6) 
-    
-    # prevalence
-    I_NU = 1-exp(-rNU)
-    I_ND = 1-exp(-rND)
-    I_SU = 1-exp(-rSU)
-    I_SD = 1-exp(-rSD)
-    I_AU = 1-exp(-rAU)
-    I_AD = 1-exp(-rAD)
-    
-    # proportion of hosts that have no parasites from W or 
-    p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) 
-    p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) 
-    p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) 
-    p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) 
-    p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) 
-    p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6))) 
-    
-    # contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite
-    ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2]))
-    ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2]))
-    ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2]))
-    ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2]))
-    ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2]))
-    ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2]))
-    
-    # converted contribution of W or R with mixed infection cos----
-    h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
-    h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost))
-    h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
-    h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
-    h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
-    h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost))
-    
-    #transmissibility is reduced to the number of infected hosts, 
-    #receiving hosts needs to have open "resource room" for accepting new strains
-    B_PW = b*(
-      h_NU[1] +
-        h_ND[1] + 
-        h_SU[1] +
-        h_SD[1] +
-        h_AU[1] +
-        h_AD[1] 
-    )
-    
-    
-    B_PR = b*(
-      h_NU[2]+
-        h_ND[2]+ 
-        h_SU[2]+
-        h_SD[2]+
-        h_AU[2]+
-        h_AD[2]
-    )
-    
-    fracN = P_acc_N/(NU+ND+1e-6)
-    fracS = P_acc_S/(SU+SD+1e-6)
-    fracA = P_acc_A/(AU+AD+1e-6)
-    
-    # per capita biting rate for hosts who are not treated
-    Bu = (B_PW+B_PR)
-    # per capita biting rate for hosts who are drug treated
-    Bd = B_PR
-    #print(B)
-    # probability that the infected strain has not been seen by the adaptive immune system
-    #infN = (1-1/nstrains)^fracN
-    infN = 1 #for generalized immunity model, infectivity values are fixed over time
-    #infS = (1-1/nstrains)^fracS
-    infS = 1#red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5
-    #infA = (1-1/nstrains)^fracA
-    infA = 1
-    #print(infN)
-    #naive class getting infected, and treated
-    NU2ND = (Bu*(1-rNU/capacity)*infN*d1)
-    
-    #naive U class getting enough # of infections to move to SU
-    NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN))))
-    #naive U class die from infections
-    NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU_gen_N*muNU)
-    
-    #naive D class getting enough # of infections to move to SD
-    ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN))))
-    
-    #naive D class die from infections
-    NDInfDie = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND)
-    
-    #SU class getting infected, show symptoms, and treated
-    SU2SD = (Bu*(1-rSU/capacity)*infS*d2)
-    
-    #SU class getting enough # of infections to move to AU
-    SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS))))
-    
-    #SD class getting enough # of infections to move to AD
-    SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS))))
-    
-    #AU class getting infected, show symptoms, and treated
-    AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)
-    
-    muWU_gen_N = muWU*exp(-red*P_acc_N)+1
-    muWD_gen_N = muWD*exp(-red*P_acc_N)+1
-    muRU_gen_N = muRU*exp(-red*P_acc_N)+1
-    muRD_gen_N = muRD*exp(-red*P_acc_N)+1
-    
-    muWU_gen_S = muWU*exp(-red*P_acc_S)+1
-    muWD_gen_S = muWD*exp(-red*P_acc_S)+1
-    muRU_gen_S = muRU*exp(-red*P_acc_S)+1
-    muRD_gen_S = muRD*exp(-red*P_acc_S)+1
-    
-    muWU_gen_A = muWU*exp(-red*P_acc_A)+1
-    muWD_gen_A = muWD*exp(-red*P_acc_A)+1
-    muRU_gen_A = muRU*exp(-red*P_acc_A)+1
-    muRD_gen_A = muRD*exp(-red*P_acc_A)+1
-    
-    # Parasite pop ODE------
-    ## sensitive-------
-    dPW_NU.dt = 0
-    
-    #print(dNU.dt)
-    
-    dPW_ND.dt = 0
-    
-    #print(dND.dt)
-    
-    dPW_SU.dt = 0
-    
-    dPW_SD.dt = 0 
-    
-    dPW_AU.dt = 0
-    
-    dPW_AD.dt = 0 
-    
-    #resistant--
-    dPR_NU.dt = B_PR*NU*infN*(1-rNU/capacity) + #new infections
-      (PR_ND/Tau)- # drug loose effectiveness
-      (PR_NU*NU2ND) -  #naive class getting infected, and treated
-      (PR_NU*NU2SU) - #naive class getting enough # of infections to move to S
-      (PR_NU*NUInfDie)- #naive class die from infections
-      (PR_NU*att)- #normal death rate
-      (PR_NU*muRU_gen_N)#clearance of infection
-    
-    #print(dNU.dt)
-    
-    dPR_ND.dt = B_PR*ND*infN*(1-rND/capacity) + #new infections
-      (PR_NU*NU2ND) - #naive class getting infected, and treated
-      (PR_ND/Tau) - # drug loose effectiveness
-      (PR_ND*ND2SD) -#naive class getting enough # of infections to move to S
-      (PR_ND*NDInfDie) -#naive class die from infections
-      (PR_ND*att)- #normal death rate
-      (PR_ND*muRD_gen_N)#clearance from drug
-    
-    #print(dND.dt)
-    
-    dPR_SU.dt = B_PR*SU*infS*(1-rSU/capacity) + #new infections
-      (PR_NU*NU2SU) +#naive class getting enough # of infections to move to S
-      (PR_SD/Tau)-  # drug loose effectiveness
-      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PR_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
-      (PR_SU*att)- #normal death rate
-      (PR_SU*muRU_gen_S)#clearance of infection
-    
-    dPR_SD.dt = B_PR*SD*infS*(1-rSD/capacity) + #new infections
-      (PR_ND*ND2SD)+ #naive class getting enough # of infections to move to S
-      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PR_SD/Tau) - # drug loose effectiveness
-      (PR_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
-      (PR_SD*att)- #normal death rate
-      (PR_SD*muRD_gen_S)#clearance from drug
-    
-    dPR_AU.dt = B_PR*AU*infA*(1-rAU/capacity) + #new infections
-      (PR_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
-      (PR_AD/Tau)-  # drug loose effectiveness
-      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PR_AU*att)- #normal death rate
-      (PR_AU*muRU_gen_A)#clearance from infection
-    
-    dPR_AD.dt = B_PR*AD*infA*(1-rAD/capacity) + #new infections
-      (PR_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
-      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PR_AD/Tau)-  # drug loose effectiveness
-      (PR_AD*att)- #normal death rate
-      (PR_AD*muRD_gen_A)#clearance from drug
-    
-    # Generalized immunity ODE------
-    dNU.dt = birth + 
-      (ND/Tau)- # drug loose effectiveness
-      NU*NU2ND -  #naive class getting infected, and treated
-      NU*NU2SU - #naive class getting enough # of infections to move to S
-      NU*NUInfDie- #naive class die from infections
-      (NU*att)#normal death rate
-    
-    
-    #print(dNU.dt)
-    
-    dND.dt = NU*NU2ND - #naive class getting infected, and treated
-      (ND/Tau) - # drug loose effectiveness
-      ND*ND2SD -#naive class getting enough # of infections to move to S
-      ND*NDInfDie -#naive class die from infections
-      (ND*att) #normal death rate
-    #print(dND.dt)
-    
-    dSU.dt =  NU*NU2SU +#naive class getting enough # of infections to move to S
-      (SD/Tau)-  # drug loose effectiveness
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      SU*SU2AU -#symptomatic class getting enough # of infections to move to AU
-      (SU*att) #normal death rate
-    
-    dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      (SD/Tau) - # drug loose effectiveness
-      SD*SD2AD- #symptomatic class getting enough # of infections to move to AD
-      (SD*att) #normal death rate
-    
-    dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS
-      (AD/Tau)-  # drug loose effectiveness
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AU*att) #normal death rate
-    
-    dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AD/Tau)-  # drug loose effectiveness
-      (AD*att) #normal death rate
-    
-    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
-    
-    
-    # track the mean number of times infected per immunity class------
-    dP_acc_N.dt = NU*Bu*(1-rNU/capacity)*(1-infN) + ND*Bd*(1-rND/capacity)*(1-infN)+
-      NU*Bu*(1-rNU/capacity)*infN*muWU_gen_N + ND*Bd*(1-rND/capacity)*infN*muRD_gen_N- # total number of infections received
-      (NU*(att+NUInfDie))*fracN- #when naive U class die how many infections deducted
-      (ND*(att+NDInfDie))*fracN- #when naive D class die many infections deducted
-      (NU*NU2SU+ND*ND2SD)*fracN #when naive class leave to S
-    
-    dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S
-      SU*Bu*(1-rSU/capacity)*(infS*muWU_gen_S+(1-infS)) + 
-      SD*Bd*(1-rSD/capacity)*(infS*muRD_gen_S+(1-infS))- #total number of new infections
-      att*P_acc_S- # normal death
-      (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic
-    
-    dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic
-      AU*Bu*(1-rAU/capacity)*(infA*muWU_gen_A+(1-infA)) + 
-      AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections
-      (ImmLost+att)*P_acc_A # normal death
-    
-    
-    a <- 365/(2*pi)
-    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
-    #db.dt =0
-    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
-    
-    return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
-                  dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
-                  dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
-                  dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt)))
-    })
-}
-
-SEAI_wildOnly<-function(t,y,p){
-  
-  with(as.list(c(y,p)),{
-    # total host population size
-    N = NU+ND+SU+SD+AU+AD
-    # total parasite population size
-    #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD
-    #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD
-    #P = PW+PR
-    # parasite/host ratio
-    rNU <-  (PW_NU+PR_NU)/(NU+1e-6) 
-    rND <-  (PW_ND+PR_ND)/(ND+1e-6) 
-    rSU <-  (PW_SU+PR_SU)/(SU+1e-6) 
-    rSD <-  (PW_SD+PR_SD)/(SD+1e-6) 
-    rAU <-  (PW_AU+PR_AU)/(AU+1e-6) 
-    rAD <-  (PW_AD+PW_AD)/(AD+1e-6) 
-    
-    # prevalence
-    I_NU = 1-exp(-rNU)
-    I_ND = 1-exp(-rND)
-    I_SU = 1-exp(-rSU)
-    I_SD = 1-exp(-rSD)
-    I_AU = 1-exp(-rAU)
-    I_AD = 1-exp(-rAD)
-    
-    # proportion of hosts that have no parasites from W or 
-    p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) 
-    p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) 
-    p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) 
-    p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) 
-    p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) 
-    p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6))) 
-    
-    # contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite
-    ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2]))
-    ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2]))
-    ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2]))
-    ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2]))
-    ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2]))
-    ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2]))
-    
-    # converted contribution of W or R with mixed infection cos----
-    h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
-    h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost))
-    h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
-    h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
-    h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
-    h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost))
-    
-    #transmissibility is reduced to the number of infected hosts, 
-    #receiving hosts needs to have open "resource room" for accepting new strains
-    B_PW = b*(
-      h_NU[1] +
-        h_ND[1] + 
-        h_SU[1] +
-        h_SD[1] +
-        h_AU[1] +
-        h_AD[1] 
-    )
-    
-    
-    B_PR = b*(
-      h_NU[2]+
-        h_ND[2]+ 
-        h_SU[2]+
-        h_SD[2]+
-        h_AU[2]+
-        h_AD[2]
-    )
-    
-    
-    fracN = P_acc_N/(NU+ND+1e-6)
-    fracS = P_acc_S/(SU+SD+1e-6)
-    fracA = P_acc_A/(AU+AD+1e-6)
-    
-    # per capita biting rate for hosts who are not treated
-    Bu = (B_PW+B_PR)
-    # per capita biting rate for hosts who are drug treated
-    Bd = B_PR
-    #print(B)
-    
-    
-    # probability that the infected strain has not been seen by the adaptive immune system
-    # probability that the infected strain has not been seen by the adaptive immune system
-    # probability that the infected strain has not been seen by the adaptive immune system
-    #infN = (1-1/nstrains)^fracN
-    infN = 1 #for generalized immunity model, infectivity values are fixed over time
-    #infS = (1-1/nstrains)^fracS
-    infS = 1#red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5
-    #infA = (1-1/nstrains)^fracA
-    infA = 1
-    #print(infN)
-    #naive class getting infected, and treated
-    NU2ND = (Bu*(1-rNU/capacity)*infN*d1)
-    
-    #naive U class getting enough # of infections to move to SU
-    NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN))))
-    #naive U class die from infections
-    NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU_gen_N*muNU)
-    
-    #naive D class getting enough # of infections to move to SD
-    ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN))))
-    
-    #naive D class die from infections
-    NDInfDie = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND)
-    
-    #SU class getting infected, show symptoms, and treated
-    SU2SD = (Bu*(1-rSU/capacity)*infS*d2)
-    
-    #SU class getting enough # of infections to move to AU
-    SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS))))
-    
-    #SD class getting enough # of infections to move to AD
-    SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS))))
-    
-    #AU class getting infected, show symptoms, and treated
-    AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)
-    
-    muWU_gen_N = muWU*exp(-red*P_acc_N)+1
-    muWD_gen_N = muWD*exp(-red*P_acc_N)+1
-    muRU_gen_N = muRU*exp(-red*P_acc_N)+1
-    muRD_gen_N = muRD*exp(-red*P_acc_N)+1
-    
-    muWU_gen_S = muWU*exp(-red*P_acc_S)+1
-    muWD_gen_S = muWD*exp(-red*P_acc_S)+1
-    muRU_gen_S = muRU*exp(-red*P_acc_S)+1
-    muRD_gen_S = muRD*exp(-red*P_acc_S)+1
-    
-    muWU_gen_A = muWU*exp(-red*P_acc_A)+1
-    muWD_gen_A = muWD*exp(-red*P_acc_A)+1
-    muRU_gen_A = muRU*exp(-red*P_acc_A)+1
-    muRD_gen_A = muRD*exp(-red*P_acc_A)+1
-    
-    # Parasite pop ODE------
-    ## sensitive-------
-    dPW_NU.dt = B_PW*NU*infN*(1-rNU/capacity) + #new infections
-      (PW_ND/Tau)- # drug loose effectiveness
-      (PW_NU*NU2ND) -  #naive class getting infected, and treated
-      (PW_NU*NU2SU) - #naive class getting enough # of infections to move to S
-      (PW_NU*NUInfDie)- #naive class die from infections
-      (PW_NU*att)- #normal death rate
-      (PW_NU*muWU_gen_N)#clearance of infection
-    
-    
-    #print(dNU.dt)
-    
-    dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) + #new infections
-      (PW_NU*NU2ND) - #naive class getting infected, and treated
-      (PW_ND/Tau) - # drug loose effectiveness
-      (PW_ND*ND2SD) -#naive class getting enough # of infections to move to S
-      (PW_ND*NDInfDie) -#naive class die from infections
-      (PW_ND*att)- #normal death rate
-      (PW_ND*muWD_gen_N)#clearance from drug
-    
-    #print(dND.dt)
-    
-    dPW_SU.dt = B_PW*SU*infS*(1-rSU/capacity) + #new infections
-      (PW_NU*NU2SU) +#naive class getting enough # of infections to move to S
-      (PW_SD/Tau)-  # drug loose effectiveness
-      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PW_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
-      (PW_SU*att)- #normal death rate
-      (PW_SU*muWU_gen_S)#clearance of infection
-    
-    dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections
-      (PW_ND*ND2SD)+ #naive class getting enough # of infections to move to S
-      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PW_SD/Tau) - # drug loose effectiveness
-      (PW_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
-      (PW_SD*att)- #normal death rate
-      (PW_SD*muWD_gen_S)#clearance from drug
-    
-    dPW_AU.dt = B_PW*AU*infA*(1-rAU/capacity) + #new infections
-      (PW_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
-      (PW_AD/Tau)-  # drug loose effectiveness
-      (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PW_AU*att)- #normal death rate
-      (PW_AU*muWU_gen_A)#clearance from infection
-    
-    dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) +  #new infections
-      (PW_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
-      (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PW_AD/Tau)-  # drug loose effectiveness
-      (PW_AD*att)- #normal death rate
-      (PW_AD*muWD_gen_A)#clearance from drug
-    
-    ## resistant-------
-    dPR_NU.dt = 0
-    
-    #print(dNU.dt)
-    
-    dPR_ND.dt = 0
-    
-    #print(dND.dt)
-    
-    dPR_SU.dt = 0
-    
-    dPR_SD.dt = 0
-    
-    dPR_AU.dt = 0
-    
-    dPR_AD.dt = 0
-    # Generalized immunity ODE------
-    dNU.dt = birth + 
-      (ND/Tau)- # drug loose effectiveness
-      NU*NU2ND -  #naive class getting infected, and treated
-      NU*NU2SU - #naive class getting enough # of infections to move to S
-      NU*NUInfDie- #naive class die from infections
-      (NU*att)#normal death rate
-    
-    
-    #print(dNU.dt)
-    
-    dND.dt = NU*NU2ND - #naive class getting infected, and treated
-      (ND/Tau) - # drug loose effectiveness
-      ND*ND2SD -#naive class getting enough # of infections to move to S
-      ND*NDInfDie -#naive class die from infections
-      (ND*att) #normal death rate
-    #print(dND.dt)
-    
-    dSU.dt =  NU*NU2SU +#naive class getting enough # of infections to move to S
-      (SD/Tau)-  # drug loose effectiveness
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      SU*SU2AU -#symptomatic class getting enough # of infections to move to AU
-      (SU*att) #normal death rate
-    
-    dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      (SD/Tau) - # drug loose effectiveness
-      SD*SD2AD- #symptomatic class getting enough # of infections to move to AD
-      (SD*att) #normal death rate
-    
-    dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS
-      (AD/Tau)-  # drug loose effectiveness
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AU*att) #normal death rate
-    
-    dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AD/Tau)-  # drug loose effectiveness
-      (AD*att) #normal death rate
-    
-    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
-    
-    
-    # track the mean number of times infected per immunity class------
-    dP_acc_N.dt = NU*Bu*(1-rNU/capacity)*(1-infN) + ND*Bd*(1-rND/capacity)*(1-infN)+
-      NU*Bu*(1-rNU/capacity)*infN*muWU_gen_N + ND*Bd*(1-rND/capacity)*infN*muRD_gen_N- # total number of infections received
-      (NU*(att+NUInfDie))*fracN- #when naive U class die how many infections deducted
-      (ND*(att+NDInfDie))*fracN- #when naive D class die many infections deducted
-      (NU*NU2SU+ND*ND2SD)*fracN #when naive class leave to S
-    
-    dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S
-      SU*Bu*(1-rSU/capacity)*(infS*muWU_gen_S+(1-infS)) + 
-      SD*Bd*(1-rSD/capacity)*(infS*muRD_gen_S+(1-infS))- #total number of new infections
-      att*P_acc_S- # normal death
-      (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic
-    
-    dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic
-      AU*Bu*(1-rAU/capacity)*(infA*muWU_gen_A+(1-infA)) + 
-      AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections
-      (ImmLost+att)*P_acc_A # normal death
-    
-    
-    a <- 365/(2*pi)
-    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
-    #db.dt =0
-    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
-    
-    return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
-                  dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
-                  dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
-                  dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt)))
-  })
-}
-
-allPara<-read.csv(paramFile,sep=",",header=T)
-
-currentPara<-allPara[allPara$No==No,]
-
-#t2 = seq(from=0,to=30*365,by=1) 
-#N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
-#       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
-#       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
-#       b=allPara[1,]$bavgnull*(1-allPara[1,]$seasonality))
-#for (i in 1:42){
-#  out9 = ode(y=N0,times=t2,func=SEAI_p,parms=allPara[i,])
-#  out9_freqRes[i] = (out9[,8]+out9[,9]+out9[,10]+out9[,11]+out9[,12]+out9[,13])/(out9[,2]+out9[,3]+out9[,4]+out9[,5]+out9[,6]+out9[,7]+out9[,8]+out9[,9]+out9[,10]+out9[,11]+out9[,12]+out9[,13]) 
-#}
-
-
-
-
-
-# Step 1 get no drug treating intensity, wildOnly-----
-
-currentPara$d1<-0
-currentPara$d2<-0
-currentPara$d3<-0
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
-            PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
-            NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
-            P_acc_N=0,P_acc_S=0,P_acc_A=0,
-            b=currentPara$bavgnull*(1-currentPara$seasonality))
-
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>1000), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1)))
-
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
-
-storeOut<-tail(out2,1)[,2:23]
-# Step 2 treat occasionally, wildOnly-----
-
-currentPara$d1<-0.5
-currentPara$d2<-0.5
-currentPara$d3<-0.5
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
-       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
-       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
-       P_acc_N=0,P_acc_S=0,P_acc_A=0,
-       b=currentPara$bavgnull*(1-currentPara$seasonality))
-
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1)))
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
-
-storeOut<-rbind(storeOut,tail(out2,1)[,2:23])
-# Step 3 full treatment, wildOnly-----
-
-currentPara$d1<-1
-currentPara$d2<-1
-currentPara$d3<-1
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
-       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
-       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
-       P_acc_N=0,P_acc_S=0,P_acc_A=0,
-       b=currentPara$bavgnull*(1-currentPara$seasonality))
-
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1)))
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
-
-storeOut<-rbind(storeOut,tail(out2,1)[,2:23])
-
-# Step 4 resistant strains only, no treatment  ------
-currentPara$d1<-0
-currentPara$d2<-0
-currentPara$d3<-0
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
-       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
-       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
-       P_acc_N=0,P_acc_S=0,P_acc_A=0,
-       b=currentPara$bavgnull*(1-currentPara$seasonality))
-
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_resOnly,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_resOnly,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=0,currentPara,tail(out2,1)))
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
-
-# Step 5 treat occasionally, wildOnly + resistant-----
-
-currentPara$d1<-0.5
-currentPara$d2<-0.5
-currentPara$d3<-0.5
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = unlist(storeOut[2,])
-N0["PR_NU"]<-2
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=0,currentPara,tail(out2,1)))
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
-
-# Step 6 treat fully, wildOnly + resistant-----
-
-currentPara$d1<-1
-currentPara$d2<-1
-currentPara$d3<-1
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = unlist(storeOut[3,])
-N0["PR_NU"]<-2
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=0,currentPara,tail(out2,1)))
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
diff --git a/GeneralizedImmunityExp1.R b/GeneralizedImmunityExp1.R
deleted file mode 100644
index 8a99d2b..0000000
--- a/GeneralizedImmunityExp1.R
+++ /dev/null
@@ -1,1261 +0,0 @@
-require(deSolve)
-library(tidyverse)
-args <- commandArgs(trailingOnly = TRUE)
-print(args)
-paramFile<-args[1]
-No<-as.integer(args[2])
-resultsFolder <-args[3]
-jobid<-as.integer(args[4])
-#paramFile<-"ParamListGen.csv"
-#No<-3
-#resultsFolder<-"./"
-outFile<-paste(resultsFolder,"/",strsplit(paramFile,split="[.]")[[1]][1],"_results.csv",sep="")
-
-SEAI_p<-function(t,y,p){
-  
-  with(as.list(c(y,p)),{
-    # total host population size
-    #print(t)
-    N = NU+ND+SU+SD+AU+AD
-    # total parasite population size
-    #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD
-    #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD
-    #P = PW+PR
-    
-    # parasite/host ratio
-    rNU <-  (PW_NU+PR_NU)/(NU+1e-6) 
-    rND <-  (PW_ND+PR_ND)/(ND+1e-6) 
-    rSU <-  (PW_SU+PR_SU)/(SU+1e-6) 
-    rSD <-  (PW_SD+PR_SD)/(SD+1e-6) 
-    rAU <-  (PW_AU+PR_AU)/(AU+1e-6) 
-    rAD <-  (PW_AD+PW_AD)/(AD+1e-6) 
-    
-    # prevalence
-    I_NU = 1-exp(-rNU)
-    I_ND = 1-exp(-rND)
-    I_SU = 1-exp(-rSU)
-    I_SD = 1-exp(-rSD)
-    I_AU = 1-exp(-rAU)
-    I_AD = 1-exp(-rAD)
-    
-    # proportion of hosts that have no parasites from W or 
-    p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) 
-    p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) 
-    p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) 
-    p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) 
-    p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) 
-    p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6))) 
-    
-    # contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite
-    ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2]))
-    ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2]))
-    ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2]))
-    ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2]))
-    ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2]))
-    ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2]))
-    
-    # converted contribution of W or R with mixed infection cos----
-    h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
-    h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost))
-    h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
-    h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
-    h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
-    h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost))
-    
-    #transmissibility is reduced to the number of infected hosts, 
-    #receiving hosts needs to have open "resource room" for accepting new strains
-    B_PW = b*(
-      h_NU[1] +
-        h_ND[1] + 
-        h_SU[1] +
-        h_SD[1] +
-        h_AU[1] +
-        h_AD[1] 
-    )
-    
-    
-    B_PR = b*(
-      h_NU[2]+
-        h_ND[2]+ 
-        h_SU[2]+
-        h_SD[2]+
-        h_AU[2]+
-        h_AD[2]
-    )
-    
-    
-    
-    # per capita biting rate for hosts who are not treated
-    Bu = (B_PW+B_PR)
-    # per capita biting rate for hosts who are drug treated
-    Bd = B_PR
-    #print(B)
-    
-    # pop size +1 is to prevent N = 0 for undefined scenario
-    fracN = P_acc_N/(NU+ND+1e-6)
-    fracS = P_acc_S/(SU+SD+1e-6)
-    fracA = P_acc_A/(AU+AD+1e-6)
-    
-    muWU_gen_N = 1/((1/muWU)*exp(-red*fracN)+1)
-    muWD_gen_N = 1/((1/muWD)*exp(-red*fracN)+1)
-    muRU_gen_N = 1/((1/muRU)*exp(-red*fracN)+1)
-    muRD_gen_N = 1/((1/muRD)*exp(-red*fracN)+1)
-    
-    muWU_gen_S = 1/((1/muWU)*exp(-red*fracS)+1)
-    muWD_gen_S = 1/((1/muWD)*exp(-red*fracS)+1)
-    muRU_gen_S = 1/((1/muRU)*exp(-red*fracS)+1)
-    muRD_gen_S = 1/((1/muRD)*exp(-red*fracS)+1)
-    
-    muWU_gen_A = 1/((1/muWU)*exp(-red*fracA)+1)
-    muWD_gen_A = 1/((1/muWD)*exp(-red*fracA)+1)
-    muRU_gen_A = 1/((1/muRU)*exp(-red*fracA)+1)
-    muRD_gen_A = 1/((1/muRD)*exp(-red*fracA)+1)
-    
-    
-    # probability that the infected strain has not been seen by the adaptive immune system
-    #infN = (1-1/nstrains)^fracN
-    infN = 1 #for generalized immunity model, infectivity values are fixed over time
-    #infS = (1-1/nstrains)^fracS
-    infS = 1#red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5
-    #infA = (1-1/nstrains)^fracA
-    infA = 1
-    #print(infN)
-    #naive class getting infected, and treated
-    NU2ND = rNU*(1/8)*d1
-    
-    #naive U class getting enough # of infections to move to SU
-    #NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN))))
-    NU2SU = (rNU*muRU_gen_N)*rho1
-    
-    #naive U class die from infections
-    NUInfDie = rNU*muWU_gen_N*muNU
-    
-    #naive D class getting enough # of infections to move to SD
-    #ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN))))
-    ND2SD = (rND*muRD_gen_N)*rho1
-    
-    #naive D class die from infections
-    NDInfDie = rND*muRD_gen_N*muND
-    
-    #SU class getting infected, show symptoms, and treated
-    SU2SD = (rSU*(1/8)*d2)
-    
-    #SU class getting enough # of infections to move to AU
-    #SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS))))
-    SU2AU = (rSU*muWU_gen_S)*rho2
-    
-    #SD class getting enough # of infections to move to AD
-    #SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS))))
-    SD2AD = (rSD*muWD_gen_S)*rho2
-    
-    #AU class getting infected, show symptoms, and treated
-    AU2AD = (rAU*(1/8)*d3)
-    
-
-    
-    # Parasite pop ODE------
-    ## sensitive-------
-    dPW_NU.dt = B_PW*NU*infN*(1-rNU/capacity) + #new infections
-      (PW_ND/Tau)- # drug loose effectiveness
-      (PW_NU*NU2ND) -  #naive class getting infected, and treated
-      (PW_NU*NU2SU) - #naive class getting enough # of infections to move to S
-      (PW_NU*NUInfDie)- #naive class die from infections
-      (PW_NU*att)- #normal death rate
-      (PW_NU*muWU_gen_N)#clearance of infection
-    
-    
-    #print(dNU.dt)
-    
-    dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) + #new infections
-      (PW_NU*NU2ND) - #naive class getting infected, and treated
-      (PW_ND/Tau) - # drug loose effectiveness
-      (PW_ND*ND2SD) -#naive class getting enough # of infections to move to S
-      (PW_ND*NDInfDie) -#naive class die from infections
-      (PW_ND*att)- #normal death rate
-      (PW_ND*muWD_gen_N)#clearance from drug
-    
-    #print(dND.dt)
-    
-    dPW_SU.dt = B_PW*SU*infS*(1-rSU/capacity) + #new infections
-      (PW_NU*NU2SU) +#naive class getting enough # of infections to move to S
-      (PW_SD/Tau)-  # drug loose effectiveness
-      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PW_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
-      (PW_SU*att)- #normal death rate
-      (PW_SU*muWU_gen_S)#clearance of infection
-    
-    dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections
-      (PW_ND*ND2SD)+ #naive class getting enough # of infections to move to S
-      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PW_SD/Tau) - # drug loose effectiveness
-      (PW_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
-      (PW_SD*att)- #normal death rate
-      (PW_SD*muWD_gen_S)#clearance from drug
-    
-    dPW_AU.dt = B_PW*AU*infA*(1-rAU/capacity) + #new infections
-      (PW_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
-      (PW_AD/Tau)-  # drug loose effectiveness
-      (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PW_AU*att)- #normal death rate
-      (PW_AU*muWU_gen_A)#clearance from infection
-    
-    dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) +  #new infections
-      (PW_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
-      (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PW_AD/Tau)-  # drug loose effectiveness
-      (PW_AD*att)- #normal death rate
-      (PW_AD*muWD_gen_A)#clearance from drug
-    
-    ## resistant-------
-    dPR_NU.dt = B_PR*NU*infN*(1-rNU/capacity) + #new infections
-      (PR_ND/Tau)- # drug loose effectiveness
-      (PR_NU*NU2ND) -  #naive class getting infected, and treated
-      (PR_NU*NU2SU) - #naive class getting enough # of infections to move to S
-      (PR_NU*NUInfDie)- #naive class die from infections
-      (PR_NU*att)- #normal death rate
-      (PR_NU*muRU_gen_N)#clearance of infection
-    
-    #print(dNU.dt)
-    
-    dPR_ND.dt = B_PR*ND*infN*(1-rND/capacity) + #new infections
-      (PR_NU*NU2ND) - #naive class getting infected, and treated
-      (PR_ND/Tau) - # drug loose effectiveness
-      (PR_ND*ND2SD) -#naive class getting enough # of infections to move to S
-      (PR_ND*NDInfDie) -#naive class die from infections
-      (PR_ND*att)- #normal death rate
-      (PR_ND*muRD_gen_N)#clearance from drug
-    
-    #print(dND.dt)
-    
-    dPR_SU.dt = B_PR*SU*infS*(1-rSU/capacity) + #new infections
-      (PR_NU*NU2SU) +#naive class getting enough # of infections to move to S
-      (PR_SD/Tau)-  # drug loose effectiveness
-      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PR_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
-      (PR_SU*att)- #normal death rate
-      (PR_SU*muRU_gen_S)#clearance of infection
-    
-    dPR_SD.dt = B_PR*SD*infS*(1-rSD/capacity) + #new infections
-      (PR_ND*ND2SD)+ #naive class getting enough # of infections to move to S
-      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PR_SD/Tau) - # drug loose effectiveness
-      (PR_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
-      (PR_SD*att)- #normal death rate
-      (PR_SD*muRD_gen_S)#clearance from drug
-    
-    dPR_AU.dt = B_PR*AU*infA*(1-rAU/capacity) + #new infections
-      (PR_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
-      (PR_AD/Tau)-  # drug loose effectiveness
-      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PR_AU*att)- #normal death rate
-      (PR_AU*muRU_gen_A)#clearance from infection
-    
-    dPR_AD.dt = B_PR*AD*infA*(1-rAD/capacity) + #new infections
-      (PR_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
-      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PR_AD/Tau)-  # drug loose effectiveness
-      (PR_AD*att)- #normal death rate
-      (PR_AD*muRD_gen_A)#clearance from drug
-    
-    # Generalized immunity ODE------
-    dNU.dt = birth + 
-      (ND/Tau)- # drug loose effectiveness
-      NU*NU2ND -  #naive class getting infected, and treated
-      NU*NU2SU - #naive class getting enough # of infections to move to S
-      NU*NUInfDie- #naive class die from infections
-      (NU*att)#normal death rate
-    
-    
-    #print(dNU.dt)
-    
-    dND.dt = NU*NU2ND - #naive class getting infected, and treated
-      (ND/Tau) - # drug loose effectiveness
-      ND*ND2SD -#naive class getting enough # of infections to move to S
-      ND*NDInfDie -#naive class die from infections
-      (ND*att) #normal death rate
-    #print(dND.dt)
-    
-    dSU.dt =  NU*NU2SU +#naive class getting enough # of infections to move to S
-      (SD/Tau)-  # drug loose effectiveness
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      SU*SU2AU -#symptomatic class getting enough # of infections to move to AU
-      (SU*att) #normal death rate
-    
-    dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      (SD/Tau) - # drug loose effectiveness
-      SD*SD2AD- #symptomatic class getting enough # of infections to move to AD
-      (SD*att) #normal death rate
-    
-    dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS
-      (AD/Tau)-  # drug loose effectiveness
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AU*att) #normal death rate
-    
-    dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AD/Tau)-  # drug loose effectiveness
-      (AD*att) #normal death rate
-    
-    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
-    
-    
-    # track the mean number of times infected per immunity class------
-    dP_acc_N.dt = NU*NU2SU/rho1 + ND*ND2SD/rho1- # total number of infections received
-      (NU*(NUInfDie))*fracN- #when naive U class die how many infections deducted
-      (ND*(NDInfDie))*fracN- #when naive D class die many infections deducted
-      (NU*NU2SU+ND*ND2SD)*fracN- #when naive class leave to S
-      (ImmLost+att)*P_acc_N
-    
-    dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S
-      SU*SU2AU/rho2 + 
-      SD*SD2AD/rho2 - #total number of new infections
-      (ImmLost+att)*P_acc_S- # normal death
-      (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic
-    
-    dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic
-      AU*(Bu*(1-rAU/capacity)*(1-infA)+rAU*muWU) + 
-      AD*(Bd*(1-rAD/capacity)*(1-infA)+rAD*muRD)- #total number of new infections
-      (ImmLost+att)*P_acc_A # normal death
-    
-    
-    a <- 365/(2*pi)
-    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
-    #db.dt =0
-    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
-    #out1<- list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
-    #              dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
-    #              dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
-    #              dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt))
-    #print(out1)
-    
-    return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
-                  dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
-                  dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
-                  dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt)))
-  })
-}
-
-SEAI_resOnly<-function(t,y,p){
-  
-  with(as.list(c(y,p)),{
-    # total host population size
-    N = NU+ND+SU+SD+AU+AD
-    # total parasite population size
-    #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD
-    #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD
-    #P = PW+PR
-    
-    # parasite/host ratio
-    rNU <-  (PW_NU+PR_NU)/(NU+1e-6) 
-    rND <-  (PW_ND+PR_ND)/(ND+1e-6) 
-    rSU <-  (PW_SU+PR_SU)/(SU+1e-6) 
-    rSD <-  (PW_SD+PR_SD)/(SD+1e-6) 
-    rAU <-  (PW_AU+PR_AU)/(AU+1e-6) 
-    rAD <-  (PW_AD+PW_AD)/(AD+1e-6) 
-    
-    # prevalence
-    I_NU = 1-exp(-rNU)
-    I_ND = 1-exp(-rND)
-    I_SU = 1-exp(-rSU)
-    I_SD = 1-exp(-rSD)
-    I_AU = 1-exp(-rAU)
-    I_AD = 1-exp(-rAD)
-    
-    # proportion of hosts that have no parasites from W or 
-    p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) 
-    p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) 
-    p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) 
-    p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) 
-    p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) 
-    p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6))) 
-    
-    # contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite
-    ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2]))
-    ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2]))
-    ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2]))
-    ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2]))
-    ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2]))
-    ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2]))
-    
-    # converted contribution of W or R with mixed infection cos----
-    h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
-    h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost))
-    h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
-    h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
-    h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
-    h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost))
-    
-    #transmissibility is reduced to the number of infected hosts, 
-    #receiving hosts needs to have open "resource room" for accepting new strains
-    B_PW = b*(
-      h_NU[1] +
-        h_ND[1] + 
-        h_SU[1] +
-        h_SD[1] +
-        h_AU[1] +
-        h_AD[1] 
-    )
-    
-    
-    B_PR = b*(
-      h_NU[2]+
-        h_ND[2]+ 
-        h_SU[2]+
-        h_SD[2]+
-        h_AU[2]+
-        h_AD[2]
-    )
-    
-    fracN = P_acc_N/(NU+ND+1e-6)
-    fracS = P_acc_S/(SU+SD+1e-6)
-    fracA = P_acc_A/(AU+AD+1e-6)
-    
-    muWU_gen_N = 1/((1/muWU)*exp(-red*fracN)+1)
-    muWD_gen_N = 1/((1/muWD)*exp(-red*fracN)+1)
-    muRU_gen_N = 1/((1/muRU)*exp(-red*fracN)+1)
-    muRD_gen_N = 1/((1/muRD)*exp(-red*fracN)+1)
-    
-    muWU_gen_S = 1/((1/muWU)*exp(-red*fracS)+1)
-    muWD_gen_S = 1/((1/muWD)*exp(-red*fracS)+1)
-    muRU_gen_S = 1/((1/muRU)*exp(-red*fracS)+1)
-    muRD_gen_S = 1/((1/muRD)*exp(-red*fracS)+1)
-    
-    muWU_gen_A = 1/((1/muWU)*exp(-red*fracA)+1)
-    muWD_gen_A = 1/((1/muWD)*exp(-red*fracA)+1)
-    muRU_gen_A = 1/((1/muRU)*exp(-red*fracA)+1)
-    muRD_gen_A = 1/((1/muRD)*exp(-red*fracA)+1)
-    
-    # per capita biting rate for hosts who are not treated
-    Bu = (B_PW+B_PR)
-    # per capita biting rate for hosts who are drug treated
-    Bd = B_PR
-    #print(B)
-    # probability that the infected strain has not been seen by the adaptive immune system
-    
-    # probability that the infected strain has not been seen by the adaptive immune system
-    #infN = (1-1/nstrains)^fracN
-    infN = 1 #for generalized immunity model, infectivity values are fixed over time
-    #infS = (1-1/nstrains)^fracS
-    infS = 1#red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5
-    #infA = (1-1/nstrains)^fracA
-    infA = 1
-    #print(infN)
-    #naive class getting infected, and treated
-    NU2ND = rNU*(1/8)*d1
-    
-    #naive U class getting enough # of infections to move to SU
-    #NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN))))
-    NU2SU = (rNU*muRU_gen_N)*rho1
-    
-    #naive U class die from infections
-    NUInfDie = rNU*muWU_gen_N*muNU
-    
-    #naive D class getting enough # of infections to move to SD
-    #ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN))))
-    ND2SD = (rND*muRD_gen_N)*rho1
-    
-    #naive D class die from infections
-    NDInfDie = rND*muRD_gen_N*muND
-    
-    #SU class getting infected, show symptoms, and treated
-    SU2SD = (rSU*(1/8)*d2)
-    
-    #SU class getting enough # of infections to move to AU
-    #SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS))))
-    SU2AU = (rSU*muWU_gen_S)*rho2
-    
-    #SD class getting enough # of infections to move to AD
-    #SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS))))
-    SD2AD = (rSD*muWD_gen_S)*rho2
-    
-    #AU class getting infected, show symptoms, and treated
-    AU2AD = (rAU*(1/8)*d3)
-    
-    # Parasite pop ODE------
-    ## sensitive-------
-    dPW_NU.dt = 0
-    
-    #print(dNU.dt)
-    
-    dPW_ND.dt = 0
-    
-    #print(dND.dt)
-    
-    dPW_SU.dt = 0
-    
-    dPW_SD.dt = 0 
-    
-    dPW_AU.dt = 0
-    
-    dPW_AD.dt = 0 
-    
-    #resistant--
-    dPR_NU.dt = B_PR*NU*infN*(1-rNU/capacity) + #new infections
-      (PR_ND/Tau)- # drug loose effectiveness
-      (PR_NU*NU2ND) -  #naive class getting infected, and treated
-      (PR_NU*NU2SU) - #naive class getting enough # of infections to move to S
-      (PR_NU*NUInfDie)- #naive class die from infections
-      (PR_NU*att)- #normal death rate
-      (PR_NU*muRU_gen_N)#clearance of infection
-    
-    #print(dNU.dt)
-    
-    dPR_ND.dt = B_PR*ND*infN*(1-rND/capacity) + #new infections
-      (PR_NU*NU2ND) - #naive class getting infected, and treated
-      (PR_ND/Tau) - # drug loose effectiveness
-      (PR_ND*ND2SD) -#naive class getting enough # of infections to move to S
-      (PR_ND*NDInfDie) -#naive class die from infections
-      (PR_ND*att)- #normal death rate
-      (PR_ND*muRD_gen_N)#clearance from drug
-    
-    #print(dND.dt)
-    
-    dPR_SU.dt = B_PR*SU*infS*(1-rSU/capacity) + #new infections
-      (PR_NU*NU2SU) +#naive class getting enough # of infections to move to S
-      (PR_SD/Tau)-  # drug loose effectiveness
-      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PR_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
-      (PR_SU*att)- #normal death rate
-      (PR_SU*muRU_gen_S)#clearance of infection
-    
-    dPR_SD.dt = B_PR*SD*infS*(1-rSD/capacity) + #new infections
-      (PR_ND*ND2SD)+ #naive class getting enough # of infections to move to S
-      (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PR_SD/Tau) - # drug loose effectiveness
-      (PR_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
-      (PR_SD*att)- #normal death rate
-      (PR_SD*muRD_gen_S)#clearance from drug
-    
-    dPR_AU.dt = B_PR*AU*infA*(1-rAU/capacity) + #new infections
-      (PR_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
-      (PR_AD/Tau)-  # drug loose effectiveness
-      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PR_AU*att)- #normal death rate
-      (PR_AU*muRU_gen_A)#clearance from infection
-    
-    dPR_AD.dt = B_PR*AD*infA*(1-rAD/capacity) + #new infections
-      (PR_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
-      (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PR_AD/Tau)-  # drug loose effectiveness
-      (PR_AD*att)- #normal death rate
-      (PR_AD*muRD_gen_A)#clearance from drug
-    
-    # Generalized immunity ODE------
-    dNU.dt = birth + 
-      (ND/Tau)- # drug loose effectiveness
-      NU*NU2ND -  #naive class getting infected, and treated
-      NU*NU2SU - #naive class getting enough # of infections to move to S
-      NU*NUInfDie- #naive class die from infections
-      (NU*att)#normal death rate
-    
-    
-    #print(dNU.dt)
-    
-    dND.dt = NU*NU2ND - #naive class getting infected, and treated
-      (ND/Tau) - # drug loose effectiveness
-      ND*ND2SD -#naive class getting enough # of infections to move to S
-      ND*NDInfDie -#naive class die from infections
-      (ND*att) #normal death rate
-    #print(dND.dt)
-    
-    dSU.dt =  NU*NU2SU +#naive class getting enough # of infections to move to S
-      (SD/Tau)-  # drug loose effectiveness
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      SU*SU2AU -#symptomatic class getting enough # of infections to move to AU
-      (SU*att) #normal death rate
-    
-    dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      (SD/Tau) - # drug loose effectiveness
-      SD*SD2AD- #symptomatic class getting enough # of infections to move to AD
-      (SD*att) #normal death rate
-    
-    dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS
-      (AD/Tau)-  # drug loose effectiveness
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AU*att) #normal death rate
-    
-    dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AD/Tau)-  # drug loose effectiveness
-      (AD*att) #normal death rate
-    
-    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
-    
-    
-    # track the mean number of times infected per immunity class------
-    dP_acc_N.dt = NU*NU2SU/rho1 + ND*ND2SD/rho1- # total number of infections received
-      (NU*(NUInfDie))*fracN- #when naive U class die how many infections deducted
-      (ND*(NDInfDie))*fracN- #when naive D class die many infections deducted
-      (NU*NU2SU+ND*ND2SD)*fracN- #when naive class leave to S
-      (ImmLost+att)*P_acc_N
-    
-    dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S
-      SU*SU2AU/rho2 + 
-      SD*SD2AD/rho2 - #total number of new infections
-      (ImmLost+att)*P_acc_S- # normal death
-      (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic
-    
-    dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic
-      AU*(Bu*(1-rAU/capacity)*(1-infA)+rAU*muWU) + 
-      AD*(Bd*(1-rAD/capacity)*(1-infA)+rAD*muRD)- #total number of new infections
-      (ImmLost+att)*P_acc_A # normal death
-    
-    
-    
-    a <- 365/(2*pi)
-    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
-    #db.dt =0
-    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
-    
-    return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
-                  dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
-                  dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
-                  dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt)))
-    })
-}
-
-SEAI_wildOnly<-function(t,y,p){
-  
-  with(as.list(c(y,p)),{
-    # total host population size
-    N = NU+ND+SU+SD+AU+AD
-    # total parasite population size
-    #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD
-    #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD
-    #P = PW+PR
-    # parasite/host ratio
-    rNU <-  (PW_NU+PR_NU)/(NU+1e-6) 
-    rND <-  (PW_ND+PR_ND)/(ND+1e-6) 
-    rSU <-  (PW_SU+PR_SU)/(SU+1e-6) 
-    rSD <-  (PW_SD+PR_SD)/(SD+1e-6) 
-    rAU <-  (PW_AU+PR_AU)/(AU+1e-6) 
-    rAD <-  (PW_AD+PW_AD)/(AD+1e-6) 
-    
-    # prevalence
-    I_NU = 1-exp(-rNU)
-    I_ND = 1-exp(-rND)
-    I_SU = 1-exp(-rSU)
-    I_SD = 1-exp(-rSD)
-    I_AU = 1-exp(-rAU)
-    I_AD = 1-exp(-rAD)
-    
-    # proportion of hosts that have no parasites from W or 
-    p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) 
-    p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) 
-    p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) 
-    p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) 
-    p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) 
-    p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6))) 
-    
-    # contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite
-    ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2]))
-    ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2]))
-    ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2]))
-    ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2]))
-    ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2]))
-    ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2]))
-    
-    # converted contribution of W or R with mixed infection cos----
-    h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
-    h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost))
-    h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
-    h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
-    h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
-    h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost))
-    
-    #transmissibility is reduced to the number of infected hosts, 
-    #receiving hosts needs to have open "resource room" for accepting new strains
-    B_PW = b*(
-      h_NU[1] +
-        h_ND[1] + 
-        h_SU[1] +
-        h_SD[1] +
-        h_AU[1] +
-        h_AD[1] 
-    )
-    
-    
-    B_PR = b*(
-      h_NU[2]+
-        h_ND[2]+ 
-        h_SU[2]+
-        h_SD[2]+
-        h_AU[2]+
-        h_AD[2]
-    )
-    
-    
-    fracN = P_acc_N/(NU+ND+1e-6)
-    fracS = P_acc_S/(SU+SD+1e-6)
-    fracA = P_acc_A/(AU+AD+1e-6)
-    
-    muWU_gen_N = 1/((1/muWU)*exp(-red*fracN)+1)
-    muWD_gen_N = 1/((1/muWD)*exp(-red*fracN)+1)
-    muRU_gen_N = 1/((1/muRU)*exp(-red*fracN)+1)
-    muRD_gen_N = 1/((1/muRD)*exp(-red*fracN)+1)
-    
-    muWU_gen_S = 1/((1/muWU)*exp(-red*fracS)+1)
-    muWD_gen_S = 1/((1/muWD)*exp(-red*fracS)+1)
-    muRU_gen_S = 1/((1/muRU)*exp(-red*fracS)+1)
-    muRD_gen_S = 1/((1/muRD)*exp(-red*fracS)+1)
-    
-    muWU_gen_A = 1/((1/muWU)*exp(-red*fracA)+1)
-    muWD_gen_A = 1/((1/muWD)*exp(-red*fracA)+1)
-    muRU_gen_A = 1/((1/muRU)*exp(-red*fracA)+1)
-    muRD_gen_A = 1/((1/muRD)*exp(-red*fracA)+1)
-    
-    # per capita biting rate for hosts who are not treated
-    Bu = (B_PW+B_PR)
-    # per capita biting rate for hosts who are drug treated
-    Bd = B_PR
-    #print(B)
-    
-    
-    # probability that the infected strain has not been seen by the adaptive immune system
-    # probability that the infected strain has not been seen by the adaptive immune system
-    # probability that the infected strain has not been seen by the adaptive immune system
-    
-    # probability that the infected strain has not been seen by the adaptive immune system
-    #infN = (1-1/nstrains)^fracN
-    infN = 1 #for generalized immunity model, infectivity values are fixed over time
-    #infS = (1-1/nstrains)^fracS
-    infS = 1#red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5
-    #infA = (1-1/nstrains)^fracA
-    infA = 1
-    #print(infN)
-    #naive class getting infected, and treated
-    NU2ND = rNU*(1/8)*d1
-    
-    #naive U class getting enough # of infections to move to SU
-    #NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN))))
-    NU2SU = (rNU*muRU_gen_N)*rho1
-    
-    #naive U class die from infections
-    NUInfDie = rNU*muWU_gen_N*muNU
-    
-    #naive D class getting enough # of infections to move to SD
-    #ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN))))
-    ND2SD = (rND*muRD_gen_N)*rho1
-    
-    #naive D class die from infections
-    NDInfDie = rND*muRD_gen_N*muND
-    
-    #SU class getting infected, show symptoms, and treated
-    SU2SD = (rSU*(1/8)*d2)
-    
-    #SU class getting enough # of infections to move to AU
-    #SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS))))
-    SU2AU = (rSU*muWU_gen_S)*rho2
-    
-    #SD class getting enough # of infections to move to AD
-    #SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS))))
-    SD2AD = (rSD*muWD_gen_S)*rho2
-    
-    #AU class getting infected, show symptoms, and treated
-    AU2AD = (rAU*(1/8)*d3)
-    
-    # Parasite pop ODE------
-    ## sensitive-------
-    dPW_NU.dt = B_PW*NU*infN*(1-rNU/capacity) + #new infections
-      (PW_ND/Tau)- # drug loose effectiveness
-      (PW_NU*NU2ND) -  #naive class getting infected, and treated
-      (PW_NU*NU2SU) - #naive class getting enough # of infections to move to S
-      (PW_NU*NUInfDie)- #naive class die from infections
-      (PW_NU*att)- #normal death rate
-      (PW_NU*muWU_gen_N)#clearance of infection
-    
-    
-    #print(dNU.dt)
-    
-    dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) + #new infections
-      (PW_NU*NU2ND) - #naive class getting infected, and treated
-      (PW_ND/Tau) - # drug loose effectiveness
-      (PW_ND*ND2SD) -#naive class getting enough # of infections to move to S
-      (PW_ND*NDInfDie) -#naive class die from infections
-      (PW_ND*att)- #normal death rate
-      (PW_ND*muWD_gen_N)#clearance from drug
-    
-    #print(dND.dt)
-    
-    dPW_SU.dt = B_PW*SU*infS*(1-rSU/capacity) + #new infections
-      (PW_NU*NU2SU) +#naive class getting enough # of infections to move to S
-      (PW_SD/Tau)-  # drug loose effectiveness
-      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PW_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU
-      (PW_SU*att)- #normal death rate
-      (PW_SU*muWU_gen_S)#clearance of infection
-    
-    dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections
-      (PW_ND*ND2SD)+ #naive class getting enough # of infections to move to S
-      (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated
-      (PW_SD/Tau) - # drug loose effectiveness
-      (PW_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD
-      (PW_SD*att)- #normal death rate
-      (PW_SD*muWD_gen_S)#clearance from drug
-    
-    dPW_AU.dt = B_PW*AU*infA*(1-rAU/capacity) + #new infections
-      (PW_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU
-      (PW_AD/Tau)-  # drug loose effectiveness
-      (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PW_AU*att)- #normal death rate
-      (PW_AU*muWU_gen_A)#clearance from infection
-    
-    dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) +  #new infections
-      (PW_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD
-      (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated
-      (PW_AD/Tau)-  # drug loose effectiveness
-      (PW_AD*att)- #normal death rate
-      (PW_AD*muWD_gen_A)#clearance from drug
-    
-    ## resistant-------
-    dPR_NU.dt = 0
-    
-    #print(dNU.dt)
-    
-    dPR_ND.dt = 0
-    
-    #print(dND.dt)
-    
-    dPR_SU.dt = 0
-    
-    dPR_SD.dt = 0
-    
-    dPR_AU.dt = 0
-    
-    dPR_AD.dt = 0
-    # Generalized immunity ODE------
-    dNU.dt = birth + 
-      (ND/Tau)- # drug loose effectiveness
-      NU*NU2ND -  #naive class getting infected, and treated
-      NU*NU2SU - #naive class getting enough # of infections to move to S
-      NU*NUInfDie- #naive class die from infections
-      (NU*att)#normal death rate
-    
-    
-    #print(dNU.dt)
-    
-    dND.dt = NU*NU2ND - #naive class getting infected, and treated
-      (ND/Tau) - # drug loose effectiveness
-      ND*ND2SD -#naive class getting enough # of infections to move to S
-      ND*NDInfDie -#naive class die from infections
-      (ND*att) #normal death rate
-    #print(dND.dt)
-    
-    dSU.dt =  NU*NU2SU +#naive class getting enough # of infections to move to S
-      (SD/Tau)-  # drug loose effectiveness
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      SU*SU2AU -#symptomatic class getting enough # of infections to move to AU
-      (SU*att) #normal death rate
-    
-    dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S
-      SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated
-      (SD/Tau) - # drug loose effectiveness
-      SD*SD2AD- #symptomatic class getting enough # of infections to move to AD
-      (SD*att) #normal death rate
-    
-    dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS
-      (AD/Tau)-  # drug loose effectiveness
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AU*att) #normal death rate
-    
-    dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS
-      AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated
-      (AD/Tau)-  # drug loose effectiveness
-      (AD*att) #normal death rate
-    
-    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
-    
-    
-    # track the mean number of times infected per immunity class------
-    dP_acc_N.dt = NU*NU2SU/rho1 + ND*ND2SD/rho1- # total number of infections received
-      (NU*(NUInfDie))*fracN- #when naive U class die how many infections deducted
-      (ND*(NDInfDie))*fracN- #when naive D class die many infections deducted
-      (NU*NU2SU+ND*ND2SD)*fracN- #when naive class leave to S
-      (ImmLost+att)*P_acc_N
-    
-    dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S
-      SU*SU2AU/rho2 + 
-      SD*SD2AD/rho2 - #total number of new infections
-      (ImmLost+att)*P_acc_S- # normal death
-      (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic
-    
-    dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic
-      AU*(Bu*(1-rAU/capacity)*(1-infA)+rAU*muWU) + 
-      AD*(Bd*(1-rAD/capacity)*(1-infA)+rAD*muRD)- #total number of new infections
-      (ImmLost+att)*P_acc_A # normal death
-    
-    
-    a <- 365/(2*pi)
-    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
-    #db.dt =0
-    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
-    
-    return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
-                  dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
-                  dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
-                  dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt)))
-  })
-}
-
-allPara<-read.csv(paramFile,sep=",",header=T)
-
-currentPara<-allPara[allPara$No==No,]
-
-#t2 = seq(from=0,to=30*365,by=1) 
-#N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
-#       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
-#       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
-#       b=allPara[1,]$bavgnull*(1-allPara[1,]$seasonality))
-#for (i in 1:42){
-#  out9 = ode(y=N0,times=t2,func=SEAI_p,parms=allPara[i,])
-#  out9_freqRes[i] = (out9[,8]+out9[,9]+out9[,10]+out9[,11]+out9[,12]+out9[,13])/(out9[,2]+out9[,3]+out9[,4]+out9[,5]+out9[,6]+out9[,7]+out9[,8]+out9[,9]+out9[,10]+out9[,11]+out9[,12]+out9[,13]) 
-#}
-
-
-
-
-# Step 1 get no drug treating intensity, wildOnly-----
-
-currentPara$d1<-0
-currentPara$d2<-0
-currentPara$d3<-0
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
-            PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
-            NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
-            P_acc_N=0,P_acc_S=0,P_acc_A=0,
-            b=currentPara$bavgnull*(1-currentPara$seasonality))
-
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>1000), id.vars=1)
-ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1)))
-
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
-
-storeOut<-tail(out2,1)[,2:23]
-# Step 2 treat occasionally, wildOnly-----
-
-currentPara$d1<-0.5
-currentPara$d2<-0.5
-currentPara$d3<-0.5
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
-       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
-       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
-       P_acc_N=0,P_acc_S=0,P_acc_A=0,
-       b=currentPara$bavgnull*(1-currentPara$seasonality))
-
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1)))
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
-
-storeOut<-rbind(storeOut,tail(out2,1)[,2:23])
-# Step 3 full treatment, wildOnly-----
-
-currentPara$d1<-1
-currentPara$d2<-1
-currentPara$d3<-1
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
-       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
-       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
-       P_acc_N=0,P_acc_S=0,P_acc_A=0,
-       b=currentPara$bavgnull*(1-currentPara$seasonality))
-
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1)))
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
-
-storeOut<-rbind(storeOut,tail(out2,1)[,2:23])
-
-# Step 4 resistant strains only, no treatment  ------
-currentPara$d1<-0
-currentPara$d2<-0
-currentPara$d3<-0
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
-       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
-       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
-       P_acc_N=0,P_acc_S=0,P_acc_A=0,
-       b=currentPara$bavgnull*(1-currentPara$seasonality))
-
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_resOnly,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_resOnly,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=0,currentPara,tail(out2,1)))
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
-
-# Step 5 treat occasionally, wildOnly + resistant-----
-
-currentPara$d1<-0.5
-currentPara$d2<-0.5
-currentPara$d3<-0.5
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = unlist(storeOut[2,])
-N0["PR_NU"]<-2
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=0,currentPara,tail(out2,1)))
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
-
-# Step 6 treat fully, wildOnly + resistant-----
-
-currentPara$d1<-1
-currentPara$d2<-1
-currentPara$d3<-1
-
-t2 = seq(from=0,to=30*365,by=1) 
-N0 = unlist(storeOut[3,])
-N0["PR_NU"]<-2
-print(N0)
-out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
-
-out2<-as_tibble(out2)
-sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-  group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-print(tail(sumOut2))
-
-print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
-print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
-
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
-      (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
-  print(out2[30*365,2:23])
-  t2<-t2+30*365
-  N0 = as.numeric(out2[30*365+1,2:23])
-  names(N0)<-colnames(out2[,2:23])
-  out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
-  out2<-as_tibble(out2)
-  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
-    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
-}
-out2<-as.data.frame(out2)
-out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
-                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
-                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
-
-##you can check the output through this graph
-#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-#  facet_wrap(.~variable,scales="free_y")
-outData<-as.data.frame(cbind(jobid=jobid,wildOnly=0,currentPara,tail(out2,1)))
-if(file.exists(outFile)){
-  write.table(outData,
-              file = outFile, row.names=F, col.names=F,sep=",",
-              quote=F,append=T)
-}else{
-  write.csv(outData,
-            file = outFile, row.names=F, 
-            quote=F)
-  
-}
diff --git a/batchRuns/GeneralizedImmunityExp.R b/batchRuns/GeneralizedImmunityExp.R
index d7bdea5..1392ec1 100644
--- a/batchRuns/GeneralizedImmunityExp.R
+++ b/batchRuns/GeneralizedImmunityExp.R
@@ -1,13 +1,13 @@
 require(deSolve)
 library(tidyverse)
-#args <- commandArgs(trailingOnly = TRUE)
-#print(args)
-#paramFile<-args[1]
-#No<-as.integer(args[2])
-#resultsFolder <-args[3]
-#jobid<-as.integer(args[4])
-paramFile<-"ParamListGen.csv"
-No<-3
+args <- commandArgs(trailingOnly = TRUE)
+print(args)
+paramFile<-args[1]
+No<-as.integer(args[2])
+resultsFolder <-args[3]
+jobid<-as.integer(args[4])
+#paramFile<-"paramList.csv"
+#No<-1
 #resultsFolder<-"./"
 outFile<-paste(resultsFolder,"/",strsplit(paramFile,split="[.]")[[1]][1],"_results.csv",sep="")
 
@@ -96,22 +96,6 @@ SEAI_p<-function(t,y,p){
     fracS = P_acc_S/(SU+SD+1e-6)
     fracA = P_acc_A/(AU+AD+1e-6)
     
-    muWU_gen_N = 1/((1/muWU)*exp(-red*(fracN-((2/3)*muWU)))+1)
-    muWD_gen_N = 1/((1/muWD)*exp(-red*(fracN-((2/3)*muWD)))+1)
-    muRU_gen_N = 1/((1/muRU)*exp(-red*(fracN-((2/3)*muRU)))+1)
-    muRD_gen_N = 1/((1/muRD)*exp(-red*(fracN-((2/3)*muRD)))+1)
-    
-    muWU_gen_S = 1/((1/muWU)*exp(-red*(fracS-((2/3)*muWU)))+1)
-    muWD_gen_S = 1/((1/muWD)*exp(-red*(fracS-((2/3)*muWD)))+1)
-    muRU_gen_S = 1/((1/muRU)*exp(-red*(fracS-((2/3)*muRU)))+1)
-    muRD_gen_S = 1/((1/muRD)*exp(-red*(fracS-((2/3)*muRD)))+1)
-    
-    muWU_gen_A = 1/((1/muWU)*exp(-red*(fracA-((2/3)*muWU)))+1)
-    muWD_gen_A = 1/((1/muWD)*exp(-red*(fracA-((2/3)*muWD)))+1)
-    muRU_gen_A = 1/((1/muRU)*exp(-red*(fracA-((2/3)*muRU)))+1)
-    muRD_gen_A = 1/((1/muRD)*exp(-red*(fracA-((2/3)*muRD)))+1)
-    
-    
     # probability that the infected strain has not been seen by the adaptive immune system
     #infN = (1-1/nstrains)^fracN
     infN = 1 #for generalized immunity model, infectivity values are fixed over time
@@ -146,7 +130,20 @@ SEAI_p<-function(t,y,p){
     #AU class getting infected, show symptoms, and treated
     AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)
     
-
+    muWU_gen_N = muWU*exp(-red*P_acc_N)+1
+    muWD_gen_N = muWD*exp(-red*P_acc_N)+1
+    muRU_gen_N = muRU*exp(-red*P_acc_N)+1
+    muRD_gen_N = muRD*exp(-red*P_acc_N)+1
+    
+    muWU_gen_S = muWU*exp(-red*P_acc_S)+1
+    muWD_gen_S = muWD*exp(-red*P_acc_S)+1
+    muRU_gen_S = muRU*exp(-red*P_acc_S)+1
+    muRD_gen_S = muRD*exp(-red*P_acc_S)+1
+    
+    muWU_gen_A = muWU*exp(-red*P_acc_A)+1
+    muWD_gen_A = muWD*exp(-red*P_acc_A)+1
+    muRU_gen_A = muRU*exp(-red*P_acc_A)+1
+    muRD_gen_A = muRD*exp(-red*P_acc_A)+1
     
     # Parasite pop ODE------
     ## sensitive-------
@@ -404,20 +401,6 @@ SEAI_resOnly<-function(t,y,p){
     fracN = P_acc_N/(NU+ND+1e-6)
     fracS = P_acc_S/(SU+SD+1e-6)
     fracA = P_acc_A/(AU+AD+1e-6)
-    muWU_gen_N = 1/((1/muWU)*exp(-red*(fracN-((2/3)*muWU)))+1)
-    muWD_gen_N = 1/((1/muWD)*exp(-red*(fracN-((2/3)*muWD)))+1)
-    muRU_gen_N = 1/((1/muRU)*exp(-red*(fracN-((2/3)*muRU)))+1)
-    muRD_gen_N = 1/((1/muRD)*exp(-red*(fracN-((2/3)*muRD)))+1)
-    
-    muWU_gen_S = 1/((1/muWU)*exp(-red*(fracS-((2/3)*muWU)))+1)
-    muWD_gen_S = 1/((1/muWD)*exp(-red*(fracS-((2/3)*muWD)))+1)
-    muRU_gen_S = 1/((1/muRU)*exp(-red*(fracS-((2/3)*muRU)))+1)
-    muRD_gen_S = 1/((1/muRD)*exp(-red*(fracS-((2/3)*muRD)))+1)
-    
-    muWU_gen_A = 1/((1/muWU)*exp(-red*(fracA-((2/3)*muWU)))+1)
-    muWD_gen_A = 1/((1/muWD)*exp(-red*(fracA-((2/3)*muWD)))+1)
-    muRU_gen_A = 1/((1/muRU)*exp(-red*(fracA-((2/3)*muRU)))+1)
-    muRD_gen_A = 1/((1/muRD)*exp(-red*(fracA-((2/3)*muRD)))+1)
     
     # per capita biting rate for hosts who are not treated
     Bu = (B_PW+B_PR)
@@ -458,6 +441,20 @@ SEAI_resOnly<-function(t,y,p){
     #AU class getting infected, show symptoms, and treated
     AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)
     
+    muWU_gen_N = muWU*exp(-red*P_acc_N)+1
+    muWD_gen_N = muWD*exp(-red*P_acc_N)+1
+    muRU_gen_N = muRU*exp(-red*P_acc_N)+1
+    muRD_gen_N = muRD*exp(-red*P_acc_N)+1
+    
+    muWU_gen_S = muWU*exp(-red*P_acc_S)+1
+    muWD_gen_S = muWD*exp(-red*P_acc_S)+1
+    muRU_gen_S = muRU*exp(-red*P_acc_S)+1
+    muRD_gen_S = muRD*exp(-red*P_acc_S)+1
+    
+    muWU_gen_A = muWU*exp(-red*P_acc_A)+1
+    muWD_gen_A = muWD*exp(-red*P_acc_A)+1
+    muRU_gen_A = muRU*exp(-red*P_acc_A)+1
+    muRD_gen_A = muRD*exp(-red*P_acc_A)+1
     
     # Parasite pop ODE------
     ## sensitive-------
@@ -677,21 +674,6 @@ SEAI_wildOnly<-function(t,y,p){
     fracS = P_acc_S/(SU+SD+1e-6)
     fracA = P_acc_A/(AU+AD+1e-6)
     
-    muWU_gen_N = 1/((1/muWU)*exp(-red*(fracN-((2/3)*muWU)))+1)
-    muWD_gen_N = 1/((1/muWD)*exp(-red*(fracN-((2/3)*muWD)))+1)
-    muRU_gen_N = 1/((1/muRU)*exp(-red*(fracN-((2/3)*muRU)))+1)
-    muRD_gen_N = 1/((1/muRD)*exp(-red*(fracN-((2/3)*muRD)))+1)
-    
-    muWU_gen_S = 1/((1/muWU)*exp(-red*(fracS-((2/3)*muWU)))+1)
-    muWD_gen_S = 1/((1/muWD)*exp(-red*(fracS-((2/3)*muWD)))+1)
-    muRU_gen_S = 1/((1/muRU)*exp(-red*(fracS-((2/3)*muRU)))+1)
-    muRD_gen_S = 1/((1/muRD)*exp(-red*(fracS-((2/3)*muRD)))+1)
-    
-    muWU_gen_A = 1/((1/muWU)*exp(-red*(fracA-((2/3)*muWU)))+1)
-    muWD_gen_A = 1/((1/muWD)*exp(-red*(fracA-((2/3)*muWD)))+1)
-    muRU_gen_A = 1/((1/muRU)*exp(-red*(fracA-((2/3)*muRU)))+1)
-    muRD_gen_A = 1/((1/muRD)*exp(-red*(fracA-((2/3)*muRD)))+1)
-    
     # per capita biting rate for hosts who are not treated
     Bu = (B_PW+B_PR)
     # per capita biting rate for hosts who are drug treated
@@ -735,6 +717,20 @@ SEAI_wildOnly<-function(t,y,p){
     #AU class getting infected, show symptoms, and treated
     AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)
     
+    muWU_gen_N = muWU*exp(-red*P_acc_N)+1
+    muWD_gen_N = muWD*exp(-red*P_acc_N)+1
+    muRU_gen_N = muRU*exp(-red*P_acc_N)+1
+    muRD_gen_N = muRD*exp(-red*P_acc_N)+1
+    
+    muWU_gen_S = muWU*exp(-red*P_acc_S)+1
+    muWD_gen_S = muWD*exp(-red*P_acc_S)+1
+    muRU_gen_S = muRU*exp(-red*P_acc_S)+1
+    muRD_gen_S = muRD*exp(-red*P_acc_S)+1
+    
+    muWU_gen_A = muWU*exp(-red*P_acc_A)+1
+    muWD_gen_A = muWD*exp(-red*P_acc_A)+1
+    muRU_gen_A = muRU*exp(-red*P_acc_A)+1
+    muRD_gen_A = muRD*exp(-red*P_acc_A)+1
     
     # Parasite pop ODE------
     ## sensitive-------
@@ -882,7 +878,7 @@ SEAI_wildOnly<-function(t,y,p){
 allPara<-read.csv(paramFile,sep=",",header=T)
 
 currentPara<-allPara[allPara$No==No,]
-setwd("/home/chaillet/MalariaResistance/batchRuns")
+
 #t2 = seq(from=0,to=30*365,by=1) 
 #N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
 #       PR_NU=0,PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
@@ -898,7 +894,7 @@ setwd("/home/chaillet/MalariaResistance/batchRuns")
 
 
 # Step 1 get no drug treating intensity, wildOnly-----
-require(deSolve)
+
 currentPara$d1<-0
 currentPara$d2<-0
 currentPara$d3<-0
@@ -940,9 +936,9 @@ out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
                     P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
 
 ##you can check the output through this graph
-out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>1000), id.vars=1)
-ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-  facet_wrap(.~variable,scales="free_y")
+#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>1000), id.vars=1)
+#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
+#  facet_wrap(.~variable,scales="free_y")
 outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1)))
 
 if(file.exists(outFile)){
@@ -1000,9 +996,9 @@ out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
                     P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
 
 ##you can check the output through this graph
-out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
-ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
-  facet_wrap(.~variable,scales="free_y")
+#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1)
+#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
+#  facet_wrap(.~variable,scales="free_y")
 outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1)))
 if(file.exists(outFile)){
   write.table(outData,