Changed Generalized Immunity Script

HeLab · Apr 27, 2022 · 18fb624 · 18fb624
1 parent 4a68da4
commit 18fb624
Show file tree

Hide file tree

Showing 4 changed files with 1,389 additions and 103 deletions.
diff --git a/GeneralizedImmunityExp1.R b/GeneralizedImmunityExp1.R
@@ -6,9 +6,9 @@ paramFile<-args[1]
 No<-as.integer(args[2])
 resultsFolder <-args[3]
 jobid<-as.integer(args[4])
-#paramFile<-"paramList.csv"
-#No<-1
-#resultsFolder<-"./"
+paramFile<-"ParamListGen.csv"
+No<-3
+resultsFolder<-"./"
 outFile<-paste(resultsFolder,"/",strsplit(paramFile,split="[.]")[[1]][1],"_results.csv",sep="")
 
 SEAI_p<-function(t,y,p){
@@ -124,12 +124,15 @@ SEAI_p<-function(t,y,p){
     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))))
+    #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 = (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))))
+    #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 = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND)
@@ -138,10 +141,12 @@ SEAI_p<-function(t,y,p){
     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))))
+    #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 = (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 = (Bu*(1-rAU/capacity)*infA*omega*d3)
@@ -296,23 +301,24 @@ SEAI_p<-function(t,y,p){
 
 
     # 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_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*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/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)*(infA*muWU_gen_A+(1-infA)) + 
-      AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections
+      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
@@ -425,6 +431,8 @@ SEAI_resOnly<-function(t,y,p){
     # 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
@@ -437,12 +445,15 @@ SEAI_resOnly<-function(t,y,p){
     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))))
+    #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 = (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))))
+    #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 = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND)
@@ -451,14 +462,15 @@ SEAI_resOnly<-function(t,y,p){
     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))))
+    #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 = (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 = (Bu*(1-rAU/capacity)*infA*omega*d3)
-
+    AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)    
 
     # Parasite pop ODE------
     ## sensitive-------
@@ -573,24 +585,25 @@ SEAI_resOnly<-function(t,y,p){
 
 
     # 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_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*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/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)*(infA*muWU_gen_A+(1-infA)) + 
-      AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections
+      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
@@ -702,6 +715,8 @@ SEAI_wildOnly<-function(t,y,p){
 
     # 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
@@ -714,12 +729,15 @@ SEAI_wildOnly<-function(t,y,p){
     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))))
+    #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 = (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))))
+    #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 = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND)
@@ -728,14 +746,15 @@ SEAI_wildOnly<-function(t,y,p){
     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))))
+    #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 = (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 = (Bu*(1-rAU/capacity)*infA*omega*d3)
-
+    AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3)    
 
     # Parasite pop ODE------
     ## sensitive-------
@@ -850,21 +869,21 @@ SEAI_wildOnly<-function(t,y,p){
 
 
     # 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_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*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/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)*(infA*muWU_gen_A+(1-infA)) + 
-      AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections
+      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
 
 
@@ -897,7 +916,6 @@ currentPara<-allPara[allPara$No==No,]
 
 
 
-
 # Step 1 get no drug treating intensity, wildOnly-----
 
 currentPara$d1<-0
@@ -942,8 +960,8 @@ out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
 
 ##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")
+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)){