From 6f76ee793e7acbef156ee91df9a9672dbc90fc29 Mon Sep 17 00:00:00 2001 From: "Chaillet, Jack" Date: Wed, 10 Nov 2021 13:51:10 -0500 Subject: [PATCH] Recombination Script --- ODESRecombination-copy.R | 412 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 412 insertions(+) create mode 100644 ODESRecombination-copy.R diff --git a/ODESRecombination-copy.R b/ODESRecombination-copy.R new file mode 100644 index 0000000..d6d8da9 --- /dev/null +++ b/ODESRecombination-copy.R @@ -0,0 +1,412 @@ +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] + gamma = y[11] + with(as.list(p),{ + +#gammainit = vector() +#gamma2=vector() +#gammainit1= 0.3438303 +#numsteps = 10 +#for (i in 1:numsteps){ +#ratr1 = (PRU+PRD)/(PSU+PSD+PRU+PRD) +#P1 = ratr1 +#TauRN1=((P1*((1-d1)+(d1*TauRD1)))-(d1*TauRD1))/(1-d1) #calculate TauRN1 from equation A4 +#gammainit = gammainit1 +#for(k in 1:cmax){ +# for(l in 1:hmax+1){ +# if(h1[l,k]>c1[k]){ +# x1[l,k]=0 +# Tau1[l,k]=0 +# }else if(h1[l,k]==c1[k]){ +# x1[l,k]=1 +# Tau1[l,k]=1 +# }else{ +# x1[l,k]=(h1[l,k]*(1-S))/(((c1[k]-h1[l,k]))+(h1[l,k]*(1-S))) +# Tau1[l,k]>0 +# Tau1[l,k]=((x1[l,k])^2)+(2*x1[l,k]*(1-x1[l,k])*((1+gammainit)/2)^2)+(((1-x1[l,k])^2)*(1-(1#/(c1[k]-h1[l,k])))*(gammainit*gammainit)) +# } +# } +#} + + + + +#n = 2 + +#nprob = pbinom(h1[r,q],c1[q],prob=p) +#nprob = matrix(data=0,nrow=hmax+1,ncol=cmax) +#dn2 = ((1+gammainit)/2)^n +#Seq8=matrix(data=0,nrow=hmax+1,ncol=cmax) +#Seq9=matrix(data=0,nrow=hmax+1,ncol=cmax) +#for(q in 1:cmax){ +# for(r in 1:(hmax+1)){ +# if (h1[r,q]<=c1[q]){ +# Seq8[r,q]=(freqc1[q]) +# }else{ +# Seq8[r,q]=0 +# } +# if(h1[r,q]c1[k]){ + # x1[l,k]=0 + # Tau1[l,k]=0 + #}else if(h1[l,k]==c1[k]){ + # x1[l,k]=1 + #Tau1[l,k]=1 + #}else{ + # x1[l,k]=(h1[l,k]*(1-S))/(((c1[k]-h1[l,k]))+(h1[l,k]*(1-S))) + #Tau1[l,k]>0 + #Tau1[l,k]=((x1[l,k])^2)+(2*x1[l,k]*(1-x1[l,k])*((1+smallp)/2)^2)+(((1-x1[l,k])^2)*(1-(1/(c1[k]-h1[l,k])))*(smallp*smallp)) + #} + #} + #} + #return(((freqc1[1]*((freqh1[1,1]*Tau1[1,1])+(freqh1[2,1]*Tau1[2,1])))+(freqc1[2]*((freqh1[1,2]*Tau1[1,2])+(freqh1[2,2]*Tau1[2,2])+(freqh1[3,2]*Tau1[3,2])))+(freqc1[3]*((freqh1[1,3]*Tau1[1,3])+(freqh1[2,3]*Tau1[2,3])+(freqh1[3,3]*Tau1[3,3])+(freqh1[4,3]*Tau1[4,3])))+(freqc1[4]*((freqh1[1,4]*Tau1[1,4])+(freqh1[2,4]*Tau1[2,4])+(freqh1[3,4]*Tau1[3,4])+(freqh1[4,4]*Tau1[4,4])+(freqh1[5,4]*Tau1[5,4])))+(freqc1[5]*((freqh1[1,5]*Tau1[1,5])+(freqh1[2,5]*Tau1[2,5])+(freqh1[3,5]*Tau1[3,5])+(freqh1[4,5]*Tau1[4,5])+(freqh1[5,5]*Tau1[5,5])+(freqh1[6,5]*Tau1[6,5])))+(freqc1[6]*((freqh1[1,6]*Tau1[1,6])+(freqh1[2,6]*Tau1[2,6])+(freqh1[3,6]*Tau1[3,6])+(freqh1[4,6]*Tau1[4,6])+(freqh1[5,6]*Tau1[5,6])+(freqh1[6,6]*Tau1[6,6])+(freqh1[7,6]*Tau1[7,6]))+(freqc1[7]*((freqh1[1,7]*Tau1[1,7])+(freqh1[2,7]*Tau1[2,7])+(freqh1[3,7]*Tau1[3,7])+(freqh1[4,7]*Tau1[4,7])+(freqh1[5,7]*Tau1[5,7])+(freqh1[6,7]*Tau1[6,7])+(freqh1[7,7]*Tau1[7,7])+(freqh1[8,7]*Tau1[8,7]))+(freqc1[8]*((freqh1[1,8]*Tau1[1,8])+(freqh1[2,8]*Tau1[2,8])+(freqh1[3,8]*Tau1[3,8])+(freqh1[4,8]*Tau1[4,8])+(freqh1[5,8]*Tau1[5,8])+(freqh1[6,8]*Tau1[6,8])+(freqh1[7,8]*Tau1[7,8])+(freqh1[8,8]*Tau1[8,8])+(freqh1[9,8]*Tau1[9,8])))))/((freqc1[1])+(freqc1[2])+(freqc1[3])+(freqc1[4])+(freqc1[5])+(freqc1[6])+(freqc1[7])+(freqc1[8])))-TauRN1) + #} + #Recomb = uniroot(f=FUNC2,interval=c(0,1),S=S,Tau1=Tau1,x1=x1,c1=c1,h1=h1,freqc1=freqc1,freqh1=freqh1,f1=f1) #find root of equation 1 between 0 and 1 + #gammainit=(Recomb$root) + + #gammainit3 = 0.471481 + #gamma2final = gamma2[numsteps] + #PSU to PRU: sigma1*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PSU to PRD: (sigma1*b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PSD to PRU: (sigma1*b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PSD to PRD: (sigma1*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PRU to PSU: (sigma2*b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PRU to PSD: (sigma2*b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PRD to PSU: (sigma2*b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PRD to PSD: (sigma2*b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PSD to PSD: ((1-sigma1)*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PSU to PSU: ((1-sigma1)*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PRD to PRD: ((1-sigma2)*b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + #PRU to PRU: ((1-sigma2)*b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + dPSU.dt = -((1-sigma1)*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (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 = -((1-sigma1)*b*(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/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*(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) + #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) + dgamma.dt = (1/(PSU+PSD+PRU+PRD))*((1-gamma)/(rats))*(dPRD.dt+dPRU.dt) + return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,dgamma.dt))) + }) +} + +b = 1 + +nstrains = 150 + +del = 1 +att = 1/365/30#.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 = 5/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.15 +S = 1-wr +#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 = 0.87*0.3 +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) + +tmax=40000 +t2 = seq(from=0,to=tmax,by=10) + +P1 = 0.5 +TauRN1=((P1*((1-d1)+(d1*TauRD1)))-(d1*TauRD1))/(1-d1) #calculate TauRN1 from equation A4 +FUNC2<-function(p,E,S,Tau1,x1,f1,c1,h1,freqc1,freqh1){ + for(k in 1:cmax){ + for(l in 1:hmax+1){ + if(h1[l,k]>c1[k]){ + x1[l,k]=0 + Tau1[l,k]=0 + }else if(h1[l,k]==c1[k]){ + x1[l,k]=1 + Tau1[l,k]=1 + }else{ + x1[l,k]=(h1[l,k]*(1-S))/(((c1[k]-h1[l,k]))+(h1[l,k]*(1-S))) + Tau1[l,k]>0 + Tau1[l,k]=((x1[l,k])^2)+(2*x1[l,k]*(1-x1[l,k])*((1+p)/2)^2)+(((1-x1[l,k])^2)*(1-(1/(c1[k]-h1[l,k])))*(p*p)) + } + } + } + return(((freqc1[1]*((freqh1[1,1]*Tau1[1,1])+(freqh1[2,1]*Tau1[2,1])))+(freqc1[2]*((freqh1[1,2]*Tau1[1,2])+(freqh1[2,2]*Tau1[2,2])+(freqh1[3,2]*Tau1[3,2])))+(freqc1[3]*((freqh1[1,3]*Tau1[1,3])+(freqh1[2,3]*Tau1[2,3])+(freqh1[3,3]*Tau1[3,3])+(freqh1[4,3]*Tau1[4,3])))+(freqc1[4]*((freqh1[1,4]*Tau1[1,4])+(freqh1[2,4]*Tau1[2,4])+(freqh1[3,4]*Tau1[3,4])+(freqh1[4,4]*Tau1[4,4])+(freqh1[5,4]*Tau1[5,4])))+(freqc1[5]*((freqh1[1,5]*Tau1[1,5])+(freqh1[2,5]*Tau1[2,5])+(freqh1[3,5]*Tau1[3,5])+(freqh1[4,5]*Tau1[4,5])+(freqh1[5,5]*Tau1[5,5])+(freqh1[6,5]*Tau1[6,5])))+(freqc1[6]*((freqh1[1,6]*Tau1[1,6])+(freqh1[2,6]*Tau1[2,6])+(freqh1[3,6]*Tau1[3,6])+(freqh1[4,6]*Tau1[4,6])+(freqh1[5,6]*Tau1[5,6])+(freqh1[6,6]*Tau1[6,6])+(freqh1[7,6]*Tau1[7,6]))+(freqc1[7]*((freqh1[1,7]*Tau1[1,7])+(freqh1[2,7]*Tau1[2,7])+(freqh1[3,7]*Tau1[3,7])+(freqh1[4,7]*Tau1[4,7])+(freqh1[5,7]*Tau1[5,7])+(freqh1[6,7]*Tau1[6,7])+(freqh1[7,7]*Tau1[7,7])+(freqh1[8,7]*Tau1[8,7]))+(freqc1[8]*((freqh1[1,8]*Tau1[1,8])+(freqh1[2,8]*Tau1[2,8])+(freqh1[3,8]*Tau1[3,8])+(freqh1[4,8]*Tau1[4,8])+(freqh1[5,8]*Tau1[5,8])+(freqh1[6,8]*Tau1[6,8])+(freqh1[7,8]*Tau1[7,8])+(freqh1[8,8]*Tau1[8,8])+(freqh1[9,8]*Tau1[9,8])))))/((freqc1[1])+(freqc1[2])+(freqc1[3])+(freqc1[4])+(freqc1[5])+(freqc1[6])+(freqc1[7])+(freqc1[8])))-TauRN1) +} +Recomb = uniroot(f=FUNC2,interval=c(0,1),S=S,Tau1=Tau1,x1=x1,c1=c1,h1=h1,freqc1=freqc1,freqh1=freqh1,f1=f1) #find root of equation 1 between 0 and 1 +gammainit2=(Recomb$root) + +print(gammainit2) + +N0 = c(500,500,500,500,333,333,333,333,333,333,0.34) +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[,9])),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 = 8) +out3eq2 = vector(mode = "numeric", length = 17) +#out3eq3 = vector(mode = "numeric", length = 8) +#out3eq4 = vector(mode = "numeric", length = 8) +#out3eq5 = vector(mode="numeric", length=8) +#out3eq6 = vector(mode="numeric", length=8) +#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(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.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:17){ + #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,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])/(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]) + #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]) +} + +#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=11,ncol=17) +#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:11){ + for (j in 1:17){ + #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) \ No newline at end of file