From a619969b841a093fc03d99f9ba661c8190bab0ea Mon Sep 17 00:00:00 2001 From: Jack Chaillet Date: Fri, 3 Dec 2021 16:25:31 -0500 Subject: [PATCH] changed initial condition of sinusoidal function --- ODESRecombination-copy.R | 91 +++++++++++++++++++++++++--------------- 1 file changed, 58 insertions(+), 33 deletions(-) diff --git a/ODESRecombination-copy.R b/ODESRecombination-copy.R index c6121cf..a210dae 100644 --- a/ODESRecombination-copy.R +++ b/ODESRecombination-copy.R @@ -155,6 +155,7 @@ SEAI<-function(t,y,p){ ratr = ((PRU+PRD)*wr)/(PSU+PSD+((PRU+PRD)*wr)) ratr1 = (PRU+PRD)/(PSU+PSD+PRU+PRD) rats = (PSU+PSD)/(PSU+PSD+((PRU+PRD)*wr)) + rats1 = (PSU+PSD)/(PSU+PSD+PRU+PRD) #P1 = ratr1 #ps = 0.34 #gammainit = (ratr1*1)+(rats*ps) @@ -215,13 +216,23 @@ SEAI<-function(t,y,p){ #dPSD.dt = -((1-sigma1)*b*((PSU+PSD)/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (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 = -((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) - dPSU.dt = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSU*PSU - dPSD.dt = ((ND+SD+AD)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSD*PSD - dPRU.dt = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x^n))-muRU*PRU - dPRD.dt = ((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x^n))-muRD*PRD - e = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSD*PSD+((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x^n))-muRD*PRD - er = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x^n))-muRD*PRD - es = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*((1-x)^n))-muSD*PSD + gammavec = c((n/n-1),x) + gammavec2 = c(0,x) + if (x<0){ + gamma = max(gammavec2) + }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 + 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 #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) @@ -234,12 +245,13 @@ SEAI<-function(t,y,p){ #dgamma.dt = (1/(PSU+PSD+((1-S)*(PRU+PRD))))*((1/rats)-(gamma))*(dPRD.dt+dPRU.dt) #dgamma.dt = ratr + rats*((gamma)-ratr)-((1/(PSU+PSD+((1-S)*(PRU+PRD))))*(dPRD.dt+dPRU.dt)*((gamma)-ratr)) #dgamma.dt = -(dPRD.dt+dPRU.dt)/(dPSD.dt+dPSU.dt) - dx.dt = (1/((PSU+PSD+PRU+PRD)+(e)))*((er)-(x*(e))+((es)*((x*(PSD+PSD+PRU+PRD)-(PRU+PRD))/(PSU+PSD)))) + #dx.dt = (1/((PSU+PSD+PRU+PRD)+(e)))*((er)-(x*(e))+((es)*((x*(PSD+PSD+PRU+PRD)-(PRU+PRD))/(PSU+PSD)))) + dx.dt = (delP*(x1-x)+delR*(x-1))/((PSU+PSD)+delS) return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,dx.dt))) }) } -b = 10^(-1.5) +b = 10 nstrains = 150 @@ -258,7 +270,7 @@ omega2 = omega1 #prob of appearance of new antigen only in resistant strains L = 10000 ws = 1 -wr = 1-0.6 +wr = 1-0.1 S = 1-wr #muSU = 1 muSU=1/250#0.001 @@ -338,75 +350,88 @@ plot(out2[,1],(out2[,8]),ylim=c(0,max(out2[,9])),type="l", xlab="Time",ylab="Pop 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") +plot(out2[,1],out2[,10],ylim=c(0,max(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") +plot(out2[,1],(out2[,2]+out2[,3]+out2[,4]+out2[,5])/(out2[,6]+out2[,7]+out2[,8]+out2[,9]+out2[,10]+out2[,11])) + +M = (out2[,2]+out2[,3]+out2[,4]+out2[,5])/(out2[,6]+out2[,7]+out2[,8]+out2[,9]+out2[,10]+out2[,11]) +plot(out2[,1],0.5*(1-(((M*exp(-M)))/(1-exp(-M))))) + plot(out2[,1],out2[,12],type="l") +print(out2[,12][35000]) 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 = 17) -out3eq2 = vector(mode = "numeric", length = 11) -#out3eq3 = vector(mode = "numeric", length = 8) -#out3eq4 = vector(mode = "numeric", length = 8) -#out3eq5 = vector(mode="numeric", length=23) -#out3eq6 = vector(mode="numeric", length=8) +#out3eq = vector(mode = "numeric", length = 18) +out3eq2 = vector(mode = "numeric", length = 16) +#out3eq3 = vector(mode = "numeric", length = 18) +#out3eq4 = vector(mode = "numeric", length = 18) +#out3eq5 = vector(mode="numeric", length=18) +#out3eq7 = 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(-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)#,2,3,4) +bvalues = c(-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(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) #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) #bvalues = seq(from=-2.5,to=0.5,by=1) -wrvalues = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1) +#wrvalues = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1) +wrvalues = c(0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1) #wrvalues = c(0.20,0.22,0.24,0.26,0.28,0.30,0.32,0.34,0.36,0.38,0.40) 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.002,0.004,0.005,0.007,0.008,0.01,0.02,0.04,0.05,0.06,0.07,0.08,0.09,0.10,0.11) 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:11){ +for (i in 1:16){ #del=delvalues[i] - b=10^(bvalues[i]) + #b=10^(bvalues[i]) #wr=wrvalues[i] #d1=dvalues[i] + #d2=d1 + #d3=d1 #muSU = muSUvalues[i]*0.1 - #rho3 = rho3values[i] + rho3 = rho3values[i] + rho4=rho3 #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) out3 = ode(y=N0,times=t2,func=SEAI,parms=p2) - #out3eq[i] = out3[3999,3]+out3[3999,2] +out3[3999,4] + out3[3999,5] + #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]) #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[3999,8]+out3[3999,9])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]) - #out3eq5[i] =(out3[3999,10]+out3[3999,11])/(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]) } #plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant") #plot(bvalues,1-out4)#,ylab="Proportion Resistant") -plot(bvalues,out3eq2) +plot(rho3values,out3eq2) #print(out3eq2) print(out3eq2) -mat6 = matrix(data=0,nrow=11,ncol=11) +mat6 = matrix(data=0,nrow=16,ncol=16) #mat7 = matrix(data=0,nrow=11,ncol=17) #mat8 = matrix(data=0,nrow=11,ncol=11) #mat9 = matrix(data=0,nrow=11,ncol=18) #mat10 = matrix(data=0,nrow=11,ncol=18) #mat11 = matrix(data=0,nrow=11,ncol=18) #mat12 = matrix(data=0,nrow=11,ncol=18) -for (i in 1:11){ - for (j in 1:11){ +for (i in 1:16){ + for (j in 1:16){ #d = dvalues[i] # A = Avalues[j] - #rho3=rho3values[i] - wr = wrvalues[i] + rho3=rho3values[i] + rho4=rho3 + #wr = wrvalues[i] #muSU = muSUvalues[i]*0.01 - b = 10^(bvalues[j]) - #omega1 = omegavalues[j] + #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]