Added updated multi-locus model with seasonality

HeLab · Feb 12, 2022 · f263972 · f263972
1 parent 44ea622
commit f263972
Showing 1 changed file with 45 additions and 31 deletions.
diff --git a/ODESSeasonality (1).R b/ODESSeasonality (1).R
@@ -158,9 +158,11 @@ SEAI<-function(t,y,p){
     N = y[5]+y[6]+y[7]+y[8]+y[9]+y[10]
     m = (PSU+PSD+PRU+PRD)/N
     I = 1-exp(-m)#(1-((PSU+PSD+PRU+PRD)/(N*K))) 1-exp(-m)
-    ratr = ((PRU+PRD)*wr)/(PSU+PSD+((PRU+PRD)*wr))
+    S = 1-wr
+    h = 1-(S/(1+1*exp(-5*(I-0.5))))
+    ratr = ((PRU+PRD)*h)/(PSU+PSD+((PRU+PRD)*h))
     ratr1 = (PRU+PRD)/(PSU+PSD+PRU+PRD)
-    rats = (PSU+PSD)/(PSU+PSD+((PRU+PRD)*wr))
+    rats = (PSU+PSD)/(PSU+PSD+((PRU+PRD)*h))
     rats1 = (PSU+PSD)/(PSU+PSD+PRU+PRD)
     #P1 = ratr1
     #ps = 0.34
@@ -179,10 +181,10 @@ SEAI<-function(t,y,p){
     r = 0.5*(1-(((m*exp(-m)))/(1-exp(-m))))
     #r=0
     upsilon = (1/K)*(K-((PSU+PSD+PRU+PRD)/N))
-    delR = b*upsilon*(PRU+PRD)*wr#*I
+    delR = b*upsilon*(PRU+PRD)*h#*I
     delS = b*upsilon*(PSU+PSD)*ws#*I
     delP = delS +delR
-    n=2
+    n=3
     #TauRN1=((P1*((1-d1)+(d1*TauRD1)))-(d1*TauRD1))/(1-d1) #calculate TauRN1 from equation A4
     #FUNC2<-function(smallp,E,S,Tau1,x1,f1,c1,h1,freqc1,freqh1){
     # for(k in 1:cmax){
@@ -224,23 +226,23 @@ SEAI<-function(t,y,p){
     #dPRU.dt = -((1-sigma2)*b*((PRU+PRD)/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU)
     #dPRD.dt = -((1-sigma2)*b*((PRU+PRD)/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD)#-(del*att*(1/N)*(ND+SD+AD)*PRD)
 
-    gammavec = c((n/n-1),x)
+    gammavec = c(((n-1)/n),x)
     gammavec2 = c(0,x)
     if (x<0){
-      gamma = max(gammavec2)
+      gamma = 0
     }else{
       gamma = min(gammavec)
     }
     #x1 = ((PRU+PRD+PSU+PSD)*x-(PRU+PRD)*(1-wr))/(PSD+PSU+((PRD+PRU)*(wr)))
     #x1 = (x2^n)/(1+x2^n)
-    x1 = gamma*(rats1)+ratr1*wr
+    x1 = gamma*(rats1)+ratr1*h
     dPSU.dt = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU
     dPSD.dt = ((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD
     dPRU.dt = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU
     dPRD.dt = ((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD
-    e = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD+((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD
-    er = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD
-    es = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD
+    #e = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD+((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD
+    #er = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD
+    #es = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD
     #dNU.dt = (del*(AU+AD))+(ND/Tau)-(b*I*eps*NU*(d+rho1))-(muNU*NU) - (del*att*NU)
     dNU.dt = del+(ND/Tau)-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU)
     dND.dt = (b*I*eps*NU*d1)-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND)
@@ -275,23 +277,23 @@ rho4 = rho3#=1/nstrains #prob of immunity gain from SD to AD
 omega1 = 5/100 #prob of appearance of new antigen
 omega2 = omega1 #prob of appearance of new antigen only in resistant strains
 
-bavgnull = 0.1
-L =  0.8
+bavgnull = 10^(-1.85)
+L =  1
 kap = -1
 
 ws = 1
-wr = 1-0.03
-S = 1-wr
+wr = 1-0.6
+#S = 1-wr
 #muSU = 1
-muSU=1/250#0.001
+muSU=1/150#0.001
 #muSD = 20
 muSD=0.99
 #muRD = 1
 muRD = muSU#0.001
 #muRU = 1
 muRU = muSU#0.001
 Tau = 20#/365
-d1 = 0.87#*0.3
+d1 = 1#*0.3
 d2 = d1#/6
 d3 = d1#/10
 #m = 3
@@ -302,7 +304,7 @@ tmax=40000
 t2 = seq(from=0,to=tmax,by=1)
 
 
-N0 = c(10000,10000,1,1,333,333,333,333,333,333,0.34,bavgnull+(bavgnull*L))
+N0 = c(500,500,1,1,333,333,333,333,333,333,0.34,bavgnull+(bavgnull*L))
 out2 = ode(y=N0,times=t2,func=SEAI,parms=p2)
 
 plot(out2[,1],(out2[,2]),ylim=c(0,max(out2[,2])),type="l", xlab="Time", ylab="Population", main = paste("b=",b,", d=",d1))
@@ -327,20 +329,27 @@ lines(out2[,1],out2[,11],col="RED")
 plot(out2[,1],out2[,12],type="l")
 #print(out2[,12][35000])
 
-plot(out2[,1],abs(out2[,13]),type="l",xlim=c(0,2000))
+plot(out2[,1],abs(out2[,13]),type="l",xlim=c(39000,40000))
 
+m = (out2[,2]+out2[,3]+out2[,4]+out2[,5])/(out2[,6]+out2[,7]+out2[,8]+out2[,9]+out2[,10]+out2[,11])
 
-plot(out2[,1],(out2[,4]+out2[,5])/(out2[,2]+out2[,3]+out2[,4]+out2[,5]),type="l",xlim=c(30000,40000),ylim=c(0.75,0.85))
+plot(out2[,1],1-exp(-m),xlim=c(39000,40000))
 
-#out3eq = vector(mode = "numeric", length = 18)
-out3eq2 = vector(mode = "numeric", length = 18)
+print(mean(1-exp(-m)))
+
+plot(out2[,1],(out2[,4]+out2[,5])/(out2[,2]+out2[,3]+out2[,4]+out2[,5]),type="l",xlim=c(30000,40000))#,ylim=c(0.75,0.85))
+
+out3eq = matrix(data=0,nrow=18,ncol=366)
+#out3eq2 = vector(mode = "numeric", length = 13)
 #out3eq3 = vector(mode = "numeric", length = 18)
 #out3eq4 = vector(mode = "numeric", length = 18)
 #out3eq5 = vector(mode="numeric", length=18)
-#out3eq7 = vector(mode="numeric", length=18)
+out3eq7 = matrix(data=0,nrow=18,ncol=366)
+m = matrix(data=0,nrow=18,ncol=366)
+mavg = vector(mode="numeric",length=18)
 #out4 = vector(mode="numeric", length=22)
 #bvalues = c(0,0.0001,0.0005,0.05,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.70,0.80,0.90,1)
-bvalues = c(-3,-2.5,-2,-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.5,2)#,3,4)
+bvalues = c(-1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.5,2)#,3,4)
 bavg = vector(mode="numeric",length=13)
 #bvalues = c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2)
 kapvalues = c(-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1)
@@ -357,40 +366,45 @@ rho3values = c(0.0005,0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01
 omegavalues = c(0.01,0.012,0.014,0.016,0.018,0.02,0.022,0.024,0.026,0.028,0.030,0.032,0.034,0.036,0.038,0.040)#,0.042,0.044,0.046,0.048,0.050,0.052)
 #Lvalues = c(-3,-2.5,-2,-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.5,2)#,3,4)
 Lvalues = c(0.04,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
-bavgvalues = c(-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.5,2)
+#bavgvalues = c(-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.5,2)
+bavgvalues = c(-1.93,-1.92,-1.9,-1.8,-1.7,-1.5,-1,-0.7,-0.3,0,0.3,0.7,1,2,3,4,5,6)
+bavgvaluesnarrow = c(-1.97,-1.96,-1.95,-1.945,-1.94,-1.935,-1.93,-1.925,-1.92,-1.915,-1.91,-1.9,-1.89)
 for (i in 1:18){
   #del=delvales[i]
-  #b=10^(bvalues[i])
+  bavgnull=10^(bavgvalues[i])
   #wr=wrvalues[i]
   #d1=dvalues[i]
   #d2=d1
   #d3=d1
   #muSU = muSUvalues[i]*0.1
   #rho3 = rho3values[i]
   #rho4=rho3
-  bavgnull = 10^bavgvalues[i]
+  #bavgnull = 10^bavgvalues[i]
   #bavg[i] = bavgnull*L
   #L=Lvalues[i]
   #kap=kapvalues[i]
   #omega1 = omegavalues[i]
   p2 = list(b=b,ws=ws,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att,kap=kap,bavgnull=bavgnull)
+  N0 = c(500,500,1,1,333,333,333,333,333,333,0.34,bavgnull+(bavgnull*L))
   out3 = ode(y=N0,times=t2,func=SEAI,parms=p2)
-  #out3eq[i] = out3[39999,3]+out3[39999,2] +out3[39999,4] + out3[39999,5]
-  out3eq2[i] = (out3[39999,4]+out3[39999,5])/(out3[39999,2] + out3[39999,3] + out3[39999,4] + out3[39999,5])
+  out3eq[i,] = out3[39635:40000,3]+out3[39635:40000,2] +out3[39635:40000,4] + out3[39635:40000,5]
+  #out3eq2[i] = (out3[39999,4]+out3[39999,5])/(out3[39999,2] + out3[39999,3] + out3[39999,4] + out3[39999,5])
   #out3eq3[i] = (out3[3999,6]+out3[3999,7])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11])
   #out3eq4[i] = (out3[39999,8]+out3[39999,9])/(out3[39999,6]+out3[39999,7]+out3[39999,8]+out3[39999,9]+out3[39999,10]+out3[39999,11])
   #out3eq5[i] =(out3[39999,10]+out3[39999,11])/(out3[39999,6]+out3[39999,7]+out3[39999,8]+out3[39999,9]+out3[39999,10]+out3[39999,11])
   #out4 [i] = (1/(K*(out3[399999,6]+out3[399999,7]+out3[399999,8]+out3[399999,9]+out3[399999,10]+out3[399999,11])))*(-(((out3[3999,3]+out3[3999,2]+out3[3999,4] + out3[3999,5])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]))-K))
   #out4[i] =  (1/(K*(K)))*(-(((out3[3999,3]+out3[3999,2]+out3[3999,4] + out3[3999,5])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]))-K))
   #out3eq6[i] = (out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11])
-  #out3eq7[i] = (out3[39999,7]+out3[39999,9]+out3[39999,11])/(out3[39999,6]+out3[39999,7]+out3[39999,8]+out3[39999,9]+out3[39999,10]+out3[39999,11])
+  out3eq7[i,] = (out3[39635:40000,6]+out3[39635:40000,7]+out3[39635:40000,8]+out3[39635:40000,9]+out3[39635:40000,10]+out3[39635:40000,11])
+  m[i,] = out3eq[i,]/out3eq7[i,]
+  mavg[i] = mean(1-exp(-m[i,]))
 }
 
 #plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant")
 #plot(bvalues,1-out4)#,ylab="Proportion Resistant")
-plot(bavgvalues,out3eq2)
+plot(bavgvalues,mavg)
 #print(out3eq2)
-print(out3eq2)
+print(mavg)
 print(bavg)