From a4d4fca694fba9dbf802846e12d76d2a32058204 Mon Sep 17 00:00:00 2001 From: John Chaillet Date: Thu, 21 Oct 2021 12:47:52 -0400 Subject: [PATCH] Added constant birth rate to create stable population sizes --- ODEsImmunity.R | 69 +++++++++++++++++++++++++++++++++++++------------- 1 file changed, 51 insertions(+), 18 deletions(-) diff --git a/ODEsImmunity.R b/ODEsImmunity.R index 27568d5..ec2303e 100644 --- a/ODEsImmunity.R +++ b/ODEsImmunity.R @@ -17,29 +17,33 @@ 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) + rat = 0.8 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) dPSD.dt = -(b*(PSD/N)*ws*(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)) - (muSD)*(PSD) dPRU.dt = -(b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU) dPRD.dt = -(b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRU/N)*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD) - dNU.dt = (del*N)+ (ND/Tau)-(b*I*eps*NU*(d+rho1))-(muNU*NU) - (del*att*NU) + #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*(d+rho1))-(muNU*NU) - (del*att*NU) dND.dt = (b*I*eps*NU*d)-(ND/Tau)-(b*I*eps*ND*rho2)-(muND*ND) - (del*att*ND) dSU.dt = (b*I*eps*NU*rho1)+(SD/Tau)-(b*I*eps*SU*(d+rho3)) - (del*att*SU) dSD.dt = (b*I*eps*ND*rho2)+(b*I*eps*SU*d)-(SD/Tau)-(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*d) - (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 = (b*I*eps*SD*rho4)+(b*I*eps*omega1*AU*d)-(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 = 20000 +b = 9 nstrains = 20 -att = 0 -del = 0 -muND = 0 -muNU = 0 +att = 0.001#.000000000001 +del = 50 +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 @@ -50,8 +54,10 @@ omega2 = 1/100 #prob of appearance of new antigen only in resistant strains L = 10000 + + ws = 1 -wr = 1-0.6 +wr = 1-0.6#-0.6 #muSU = 1 muSU=0.05 #muSD = 20 @@ -61,7 +67,7 @@ muRD = 0.05 #muRU = 1 muRU = 0.05 Tau = 20#/365 -d = 0.9 +d = 0.3 #m = 3 K=10 p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att) @@ -83,36 +89,63 @@ plot(out2[,1],(out2[,4]),type="l", xlab = "Time", ylab="Population", main = "b=0 lines(out2[,1],(out2[,5]),col="RED") legend(x=40,y=6,legend=c("PRU","PRD"),bty="n",lwd = c(2,1),col=c("BLACK","RED")) -plot(out2[,1],(out2[,6]),type="l", xlab="Time", ylab="Population", main = "b=0.2239*365,d=0.30") +plot(out2[,1],(out2[,6]),ylim=c(0,max(out2[,7])),type="l", xlab="Time", ylab="Population", main = "b=0.2239*365,d=0.30") lines(out2[,1],(out2[,7]),col="RED") legend(x=40,y=1000,legend=c("NU","ND"),bty="n",lwd=c(2,1),col=c("BLACK","RED")) -plot(out2[,1],(out2[,8]),type="l", xlab="Time",ylab="Population", main = "b=0.2239*365,d=0.30") +plot(out2[,1],(out2[,8]),ylim=c(0,max(out2[,9])),type="l", xlab="Time",ylab="Population", main = "b=0.2239*365,d=0.30") lines(out2[,1],(out2[,9]),col="RED") legend(x=40,y=1000,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[,8]+out2[,9])/(out2[,10]+out2[,11])) + out3eq = vector(mode = "numeric", length = 22) out3eq2 = vector(mode = "numeric", length = 22) +out3eq3 = 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(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) -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) +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) for (i in 1:22){ - b=bvalues[i]*10 + b=bvalues[i]*1 #wr=wrvalues[i] #d=dvalues[i] - p3 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,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=p3) - #out3eq[i] = out3[3999,2] + out3[3999,3]+out3[3999,4] + out3[3999,5] + p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,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,2] + out3[3999,3]+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,10]+out3[3999,11] } -plot(bvalues,out3eq2,xlab="b/10") +plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant") +#plot(bvalues,out3eq3) print(out3eq2) +#print(out3eq3) + +mat6 = matrix(data=0,nrow=22,ncol=22) +for (i in 1:22){ + for (j in 1:22){ + #d = dvalues[i] + # A = Avalues[j] + wr = wrvalues[i] + b = bvalues[j]*10 + p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,m=m) + 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]) + + } +} +print(mat6) + +