Added most recent copy of script

HeLab · Nov 3, 2021 · 04e91be · 04e91be
1 parent 1810e5c
commit 04e91be
Showing 1 changed file with 208 additions and 0 deletions.
diff --git a/ODEsImmunity-copy.R b/ODEsImmunity-copy.R
@@ -0,0 +1,208 @@
+require(deSolve)
+#library(tidyverse)
+SEAI<-function(t,y,p){
+  PSU = y[1]
+  PSD = y[2]
+  #PS = y[3]
+  PRU = y[3]
+  PRD = y[4]
+  #PR = y[6]
+  NU = y[5]
+  ND = y[6]
+  SU = y[7]
+  SD = y[8]
+  AU = y[9]
+  AD = y[10]
+  with(as.list(p),{
+    N = y[5]+y[6]+y[7]+y[8]+y[9]+y[10]
+    m = (PSU+PSD+PRU+PRD)/N
+    I = (1-((PSU+PSD+PRU+PRD)/(N*K))) #1-exp(-m)
+    ratr = (PRU+PRD)/(PSU+PSD+PRU+PRD)
+    #rats = (PSU+PSD)/(PSU+PSD+PRU+PRD)
+    eps = (PSD+PSU+((PRD+PRU)*(wr)))/N
+    #eps = (PSD+PSU+((PRD+PRU)*(wr)))/(PSD+PSU+PRD+PRU)
+    dPSU.dt = - (b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSU)*(PSU)#-(del*att*(1/N)*(NU+SU+AU)*PSU)
+    dPSD.dt = - (b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)#-(del*att*(1/N)*(ND+SD+AD)*PSD)
+    dPRU.dt = - (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 = - (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)
+    #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)
+    dSU.dt = (b*I*eps*NU*rho1)+(SD/Tau)-(b*I*eps*SU*(d2+rho3)) - (del*att*SU)
+    dSD.dt = (ratr*b*I*eps*ND*rho2)+(b*I*eps*SU*d2)-(SD/Tau)-(ratr*b*I*eps*SD*rho4) - (del*att*SD)
+    #dAU.dt = rat*(AU+AD)*(1-((AU+AD)/L))+(AD/Tau)+(b*I*eps*SU*rho3)-(b*I*eps*omega1*AU*d) - (del*att*AU)
+    dAU.dt = (AD/Tau)+(b*I*eps*SU*rho3)-(b*I*eps*omega1*AU*d3) - (del*att*AU)
+    #dAD.dt = rat*(AU+AD)*(1-((AU+AD)/L))+(b*I*eps*SD*rho4)+(b*I*eps*omega1*AU*d)-(AD/Tau) - (del*att*AD)
+    dAD.dt = (ratr*b*I*eps*SD*rho4)+(b*I*eps*omega1*AU*d3)-(AD/Tau) - (del*att*AD)
+    return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)))
+  })
+}
+
+b = 1
+
+nstrains = 2000
+
+del = 60
+att = 1/365/del#.000000000001
+muND = 0.1
+muNU = 0.1
+rho1 = 0.5#prob of immunity gain from NU to SU
+rho2 = 0.5#prob of immunity gain from ND to SD
+rho3 = 1/nstrains #prob of immunity gain from SU to AU
+rho4 = rho3#=1/nstrains #prob of immunity gain from SD to AD
+omega1 = 1/100 #prob of appearance of new antigen
+omega2 = omega1 #prob of appearance of new antigen only in resistant strains
+
+
+L = 10000
+
+ws = 1
+wr = 1-0.6
+#muSU = 1
+muSU=1/250#0.001
+#muSD = 20
+muSD=0.99
+#muRD = 1
+muRD = muSU#0.001
+#muRU = 1
+muRU = muSU#0.001
+Tau = 20#/365
+d1 = 1
+d2 = d1#/6
+d3 = d1#/10
+#m = 3
+K=10
+p2 = list(b=b,ws=ws,wr=wr,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, K=K)
+
+
+t2 = seq(from=0,to=50000,by=10) 
+
+
+N0 = c(500,500,500,500,333,333,333,333,333,333)
+out2 = ode(y=N0,times=t2,func=SEAI,parms=p2)
+P<-rowSums(out2[,2:5])
+N<-rowSums(out2[,6:11])
+#plot(out2[,1],(out2[,2]+out2[,3]+out2[,4]+out2[,5]))
+I<-1-P/N/K
+PS<-rowSums(out2[,2:3]) 
+PR<-rowSums(out2[,4:5])
+eps<-(PS+wr*PR)/N
+plot(out2[,1],I,type="l", xlab="Time", ylab="I")
+plot(out2[,1],I*eps,type="l", xlab="Time", ylab="I*eps")
+
+plot(out2[,1],(out2[,2]),type="l", xlab="Time", ylab="Population", main = paste("b=",b,", d=",d1))
+lines(out2[,1],(out2[,3]),col="RED")
+legend("bottomright",legend=c("PSU","PSD"),bty="n",lwd = c(2,1),col=c("BLACK","RED"))
+
+plot(out2[,1],(out2[,4]),type="l", xlab = "Time", ylab="Population", main = "b=0.2239*365,d=0.30")
+lines(out2[,1],(out2[,5]),col="RED")
+legend("bottomright",legend=c("PRU","PRD"),bty="n",lwd = c(2,1),col=c("BLACK","RED"))
+
+plot(out2[,1],(out2[,12]),type="l", xlab = "Time", ylab="I", main = "b=0.2239*365,d=0.30")
+plot(out2[,1],(out2[,13]),type="l", xlab = "Time", ylab="eps", main = "b=0.2239*365,d=0.30")
+
+#parasite host ratio
+plot(out2[,1],apply(out2[,2:5],1,sum)/apply(out2[,6:11],1,sum),type="l", xlab="Time", ylab="parasite/hosts", main = paste("b=",b,", d=",d1))
+
+plot(out2[,1],(out2[,6]),ylim=c(0,max(out2[,6])),type="l", xlab="Time", ylab="Population", main = "b=0.2239*365,d=0.30")
+lines(out2[,1],(out2[,7]),col="RED")
+legend("bottomright",legend=c("NU","ND"),bty="n",lwd=c(2,1),col=c("BLACK","RED"))
+
+plot(out2[,1],(out2[,8]),ylim=c(0,max(out2[,8])),type="l", xlab="Time",ylab="Population", main = "b=0.2239*365,d=0.30")
+lines(out2[,1],(out2[,9]),col="RED")
+legend("bottomright",legend=c("SU","SD"),bty="n",lwd=c(2,1),col=c("BLACK","RED"))
+
+plot(out2[,1],out2[,10],type="l",xlab="Time",ylab="Asymptomatic Population (A)", main="b=0.2239*365,d=0.30")
+lines(out2[,1],out2[,11],col="RED")
+
+
+
+I1 =  (1/(K*(out3[,6]+out3[,7]+out3[,8]+out3[,9]+out3[,10]+out3[,11])))*(-(((out3[,3]+out3[,2]+out3[,4] + out3[,5])/(out3[,6]+out3[,7]+out3[,8]+out3[,9]+out3[,10]+out3[,11]))-K))
+
+plot(out2[,1],I1)
+
+out3eq = vector(mode = "numeric", length = 22)
+out3eq2 = vector(mode = "numeric", length = 22)
+out3eq3 = vector(mode = "numeric", length = 22)
+out3eq4 = vector(mode = "numeric", length = 22)
+out3eq5 = vector(mode="numeric", length=22)
+out3eq6 = vector(mode="numeric", length=22)
+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(-4,-3.5,-3,-2.8,-2.5,-2.2,-2,-1.9,-1.7,-1.5,-1.2,-1.1,-1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,2,4)
+#bvalues = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)
+wrvalues = c(0,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.65,0.70,0.75,0.80,0.85,0.90,1)
+dvalues = c(0,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.65,0.70,0.75,0.80,0.85,0.90,1)
+muSUvalues = c(0,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.65,0.70,0.75,0.80,0.85,0.90,1)
+delvalues = c(10,20,40,50,60,70,90,100,150,200,250,300,350,400,450,500,550,600,700,800,900,1000)
+rho3values = c(0.0005,0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.10,0.11,0.12)
+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)
+for (i in 1:22){
+  #del=delvalues[i]
+  b=10^(bvalues[i])
+  #wr=wrvalues[i]
+  #d1=dvalues[i]
+  #muSU = muSUvalues[i]*0.1
+  #rho3 = rho3values[i]
+  #omega1 = omegavalues[i]
+  p2 = list(b=b,ws=ws,wr=wr,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)
+  out3 = ode(y=N0,times=t2,func=SEAI,parms=p2)
+  out3eq[i] = out3[3999,3]+out3[3999,2] +out3[3999,4] + out3[3999,5]
+  out3eq2[i] = (out3[3999,4]+out3[3999,5])/(out3[3999,2] + out3[3999,3] + out3[3999,4] + out3[3999,5])
+  out3eq3[i] = out3[3999,6]+out3[3999,7]
+  out3eq4[i] = out3[3999,8]+out3[3999,9]
+  out3eq5[i] = out3[3999,10]+out3[3999,11]
+  out4 [i] = (1/(K*(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,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,7]+out3[3999,9]+out3[3999,11])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11])
+}
+
+#plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant")
+#plot(bvalues,1-out4)#,ylab="Proportion Resistant")
+plot(bvalues,out3eq2)
+#print(out3eq2)
+print(out3eq2)
+
+mat6 = matrix(data=0,nrow=22,ncol=22)
+#mat7 = matrix(data=0,nrow=22,ncol=22)
+#mat8 = matrix(data=0,nrow=22,ncol=22)
+#mat9 = matrix(data=0,nrow=22,ncol=22)
+#mat10 = matrix(data=0,nrow=22,ncol=22)
+#mat11 = matrix(data=0,nrow=22,ncol=22)
+#mat12 = matrix(data=0,nrow=22,ncol=22)
+for (i in 1:22){
+  for (j in 1:22){
+    #d = dvalues[i]
+    #    A = Avalues[j]
+    #rho3=rho3values[i]
+    wr = wrvalues[i]
+    #muSU = muSUvalues[i]*0.01
+    b = 10^(bvalues[j])
+    #omega1 = omegavalues[j]
+    p2 = list(b=b,ws=ws,wr=wr,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)
+    out5 = ode(y=N0,times=t2,func=SEAI,parms=p2)
+    #    PTotal = out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5]
+    #    mat[i,j] = (out5[3999,4] + out5[3999,5])/PTotal
+    #mat2[i,j] = (out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5])/6000
+    #mat3[i,j] = out5[3999,6]/out5[3999,7]
+    mat6[i,j] = (out5[3999,4]+out5[3999,5])/(out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5])
+    #mat7[i,j] = out5[3999,2]+out5[3999,3]+out5[3999,4]+out5[3999,5]
+    #mat8[i,j]=(out5[3999,10]+out5[3999,11])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
+    #mat9[i,j] = (out5[3999,8]+out5[3999,9])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
+    #mat10[i,j] = (out5[3999,6]+out5[3999,7])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
+    #mat11[i,j] = (out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
+    #mat12[i,j]=(out5[3999,7]+out5[3999,9]+out5[3999,11])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
+  }
+}
+print(mat6)
+#print(mat7)
+#print(mat8)
+#print(mat9)
+#print(mat10)
+#print(mat11)
+#print(mat12)
+
+
+
+
+