diff --git a/ODESImmunity2-copy.R b/ODESImmunity2-copy.R
new file mode 100644
index 0000000..9f35222
--- /dev/null
+++ b/ODESImmunity2-copy.R
@@ -0,0 +1,521 @@
+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]
+  beta = y[11]
+  #PMOR = y[12]
+  #HMOR = y[13]
+  #PSUACC = y[12]
+  #PSDACC = y[13]
+  #PRUACC = y[14]
+  #PRDACC = y[15]
+  #NUACC = y[18]
+  #NDACC = y[19]
+  PSUACCN2 = y[12]
+  PSDACCN2= y[13]
+  PRUACCN2 = y[14]
+  PRDACCN2 = y[15]
+  #PSUACC210 = y[16]
+  #PSDACC210 = y[17]
+  #PRUACC210 = y[18]
+  #PRDACC210 = y[19]
+  #ACC210U = y[20]
+  #ACC210D = y[21]
+  #ACC210NU = y[22]
+  #ACC210ND = y[23]
+  #PSUACCN2 = y[24]
+  #PSDACCN2 = y[25]
+  #PRUACCN2 = y[26]
+  #PRDACCN2 = y[27]
+  PSUACC = y[16]
+  PSDACC = y[17]
+  PRUACC = y[18]
+  PRDACC = y[19]
+  PSUACCA = y[20]
+  PSDACCA = y[21]
+  PRUACCA = y[22]
+  PRDACCA = y[23]
+  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-exp(-m)#(1-((PSU+PSD+PRU+PRD)/(N*K))) 1-exp(-m)
+    #pone = 1 - (m/K)
+    #ptwo = m/(K-1)
+    #pthree = m/(K-2)
+    #pfour = m/(K-3)
+    #pfive = m/(K-4)
+    #psix = m/(K-5)
+    #pseven = m/(K-6)
+    #peight = m/(K-7)
+    #pnine = m/(K-8)
+    #pten = m/(K-9)
+    ##I = pnbinom(K,size=K,prob=p) - pnbinom(0,size=K,prob=p)
+    #I = pnbinom(K,size=K,prob=pone) + pnbinom(K-1,size=K,prob=pone) + pnbinom(K-2,size=K,prob=pone) + pnbinom(K-3,size=K,prob=pone) + pnbinom(K-4,size=K,prob=pone) + pnbinom(K-5,size=K,prob=pone)+ pnbinom(K-6,size=K,prob=pone)+ pnbinom(K-7,size=K,prob=pone)+ pnbinom(K-8,size=K,prob=pone)+pnbinom(K-9,size=K,prob=pone) - pnbinom(0,size=K,prob=pone) - pnbinom(0,size=K-1,prob=ptwo) - pnbinom(0,size=K-2,prob=pthree) - pnbinom(0,size=K-3,prob=pfour) - pnbinom(0,size=K-4,prob=pfive) - pnbinom(0,size=K-5,prob=psix)- pnbinom(0,size=K-6,prob=pseven)- pnbinom(0,size=K-7,prob=peight)- pnbinom(0,size=K-8,prob=pnine)#-pnbinom(0,size=K-9,prob=pten)
+    #I = dnbinom(K,size=1,prob=pone)
+    kap = -1
+    a =365/(2*pi)
+    ps = 0
+    b = abs(beta)
+    S = 1-wr
+    k=12
+    h = 1-(S/(1+0.1*exp(-k*(I-1))))
+    # rho_inverse = vector(mode="numeric",length=nstrains)
+    # for (i in 1:nstrains){
+    #   rho_inverse[i] = 1/i
+    # }
+    #n_infections = nstrains*sum(rho_inverse)
+    #prob_infection = nstrains/n_infections
+    #prob_asymptomatic = 1/n_infections
+    #rho3=prob_asymptomatic
+    #rho4=rho3
+    #h=wr
+    #totalm = (PSUACC+PSDACC+PRUACC+PRDACC)/(K*(SU+SD))
+    frac1 = (PSUACC+PRUACC)/(SU)
+    frac2 = (PSDACC+PRDACC)/(SD)
+    frac3 = (PSUACCA+PRUACCA)/AU
+    frac4 = (PSDACCA+PRDACCA)/AD
+    #frac2 = PSDACC/(SD)
+    #frac3 = PRUACC/(SU)
+    #frac4 = PRDACC/(SD)
+    #f1 = frac1 #*(1-rho3)
+    #f2 = frac2 #*
+    #f3 = frac3
+    #f4 = frac4
+    #totalm = m
+    prob23 = min(frac1,nstrains)/nstrains
+    prob24 = min(frac2,nstrains)/nstrains
+    rhobin1 = dbinom(nstrains,size=nstrains,prob=prob23*0.999999)
+    rhobin2 = dbinom(nstrains,size=nstrains,prob=prob24*0.999999)
+    ratr = ((PRU+PRD)*h)/(PSU+PSD+((PRU+PRD)*h))
+    #rats = (PSU+PSD)/(PSU+PSD+PRU+PRD)
+    #eps = (PSD+PSU+((PRD+PRU)*(wr)))/N
+    #Ndeaths_demo = del*att*NU*t+del*att*ND*t
+    #Sdeaths = del*att*SU*t+del*att*SD*t
+    #Adeaths = (del*att*AU*t)+(del*att*AD*t)
+    #totalAhosts = AU+AD#+(del*att*AU*t)+(del*att*AD*t)
+    #rebirth = (PMOR)/(SU+SD)
+    #probbin = ((PSUACC+PSDACC+PRUACC+PRDACC))/((K)*(SU+SD))
+    #mbin = pnbinom((1/rho3),size = 1,prob=probbin,lower.tail=TRUE)
+    #meanbin = 1*(1-probbin)/probbin #mean number of treatment failures before a successful immune defense
+    #meanstrainsS = meanbin*((1-rho3)^meanbin)
+    meanstrainsS1 = (1-rho3)^(frac1)
+    meanstrainsS2 = (1-rho3)^(frac2)
+    meanstrainsS3 = (1-rho3)^(frac1)
+    meanstrainsS4 = (1-rho3)^(frac2)
+    meanstrainsS5 = (1-rho3)^(frac3)
+    meanstrainsS6 = (1-rho3)^(frac4)
+    eps = (PSD+PSU+((PRD+PRU)*(h)))/(PSD+PSU+PRD+PRU)
+    #dPSU.dt = - (b*(PSU/N)*ws*(NU+SU*(1-f1)*meanstrainsS1+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU*(1-f1)*meanstrainsS1+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSU)*(PSU)#-(del*att*(1/N)*(NU+SU+AU)*PSU)
+    #dPSU.dt = - (b*(PSU/N)*ws*(NU+SU+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSU)*(PSU)#-(del*att*(1/N)*(NU+SU+AU)*PSU)
+    dPSU.dt = - (b*(PSU/N)*ws*(NU+SU*(meanstrainsS1+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU*(meanstrainsS1+omega1)+AU*(omega1+meanstrainsS5))*(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*(1-f2)*meanstrainsS2+AD*0)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND+SD*(1-f2)*meanstrainsS2+AD*0)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)#-(del*att*(1/N)*(ND+SD+AD)*PSD)
+    #dPSD.dt = - (b*(PSU/N)*ws*(ND+SD+AD*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND+SD+AD*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)#-(del*att*(1/N)*(ND+SD+AD)*PSD)
+    dPSD.dt = - (b*(PSU/N)*ws*(ND+SD*(meanstrainsS2+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND+SD*(meanstrainsS2+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)#-(del*att*(1/N)*(ND+SD+AD)*PSD)
+    #dPRU.dt = - (b*(PRU/N)*h*(NU+SU+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU+SU+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU)
+    dPRU.dt = - (b*(PRU/N)*h*(NU+SU*(meanstrainsS3+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU+SU*(meanstrainsS3+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU)
+    #dPRD.dt = - (b*(PRU/N)*h*(ND+SD+AD*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND+SD+AD*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD)#-(del*att*(1/N)*(ND+SD+AD)*PRD)
+    dPRD.dt = - (b*(PRU/N)*h*(ND+SD*(meanstrainsS4+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND+SD*(meanstrainsS4+omega1)+AD*(omega1+meanstrainsS6))*(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*(meanstrainsS1+omega1)+(eta))) - (del*att*SU)
+    dSD.dt = (ratr*b*I*eps*ND*rho2)+(b*I*eps*SU*d2*(meanstrainsS1+omega1))-(SD/Tau)-(ratr*b*I*eps*SD*eta) - (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*eta)-(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*eta)+(b*I*eps*(omega1)*AU*d3)-(AD/Tau) - (del*att*AD)
+    dbeta.dt = ((L*bavgnull*sin(((t/a)+ps)-kap*sin((t/a)+ps))*(kap*cos((t/a)+ps)-1))/a)
+    #dPMOR.dt = (muSU)*(PSUACC)+(muSD)*PSDACC+muRU*PRUACC+muRD*PRDACC
+    #dHMOR.dt = del*att*del*t+(muNU*NUACC*I*b*eps)+(muND*ratr*NDACC*I*b*eps)
+    #dPSUACC.dt = - (b*(PSU/N)*ws*(SU*meanstrainsS1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(SU*meanstrainsS1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(b*I*eps*NU*rho1)*(PSUACCN2/NU)+(SD/Tau)*(PSDACC/SD)+(-(b*I*eps*SU*(d2+(rho3))) - (del*att*SU))*(PSUACC/SU)#*(1-rho3)
+    #dPSDACC.dt = - (b*(PSU/N)*ws*(SD*meanstrainsS2)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(SD*meanstrainsS2)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ratr*b*I*eps*ND*rho2)*(PSDACCN2/ND)+(b*I*eps*SU*d2)*(PSUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*rho4) - (del*att*SD))*(PSDACC/SD)- (muSD)*(PSDACC)#death rate included because of clearance by drugs rather than immune system
+    #dPRUACC.dt = - (b*(PRU/N)*h*(SU*meanstrainsS3)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SU*meanstrainsS3)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (b*I*eps*NU*rho1)*(PRUACCN2/NU)+(SD/Tau)*(PRDACC/SD)+(-(b*I*eps*SU*(d2+(rho3))) - (del*att*SU))*(PRUACC/SU)#*(1-rho3)
+    #dPRDACC.dt = - (b*(PRU/N)*h*(SD*meanstrainsS4)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SD*meanstrainsS4)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (ratr*b*I*eps*ND*rho2)*(PRDACCN2/ND)+(b*I*eps*SU*d2)*(PRUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*rho4) - (del*att*SD))*(PRDACC/SD)#*(1-rho4)
+    #dPSUACCN.dt = - (b*(PSU/N)*ws*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ND/Tau)*(PSDACCN/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PSUACCN/NU) - (muSU)*(PSUACCN)#*(1-rho3)
+    #dPSDACCN.dt = - (b*(PSU/N)*ws*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) +(b*I*eps*NU*d1)*(PSUACCN/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PSDACCN/ND) - (muSD)*(PSDACCN)#*(1-rho4)
+    #dPRUACCN.dt = - (b*(PRU/N)*h*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) +(ND/Tau)*(PRDACCN/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PRUACCN/NU) - (muRU)*(PRUACCN)#*(1-rho3)
+    #dPRDACCN.dt = - (b*(PRU/N)*h*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (b*I*eps*NU*d1)*(PRUACCN/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PRDACCN/ND)- (muRD)*(PRDACCN)#*(1-rho4)
+    #dPSUACC210.dt = - (b*(PSU/N)*ws*(ACC210U)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ACC210U)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(PSDACC210/Tau)-(b*I*eps*(d1))*(PSUACC210)-(muNU*PSUACCN*(ACC210NU/NU)*I*b*eps)-(muSU*PSUACC210)-(PSUACC210*(1/(8*365)))
+    #dPSDACC210.dt = - (b*(PSU/N)*ws*(ACC210D)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ACC210D)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(PSDACC210/Tau)+(b*I*eps*(d1))*(PSUACC210)-(muND*PSDACCN*(ACC210ND/ND)*I*b*eps)-(muSD*PSDACC210)-(PSDACC210*(1/(8*365)))
+    #dPRUACC210.dt = - (b*(PRU/N)*h*(ACC210U)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ACC210U)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(PRDACC210/Tau)-(b*I*eps*(d1))*(PRUACC210)-(muNU*PRUACCN*(ACC210NU/NU)*I*b*eps)-(muRU*PRUACC210)-(PRUACC210*(1/(8*365)))
+    #dPRDACC210.dt = - (b*(PRU/N)*h*(ACC210D)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ACC210D)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(PRDACC210/Tau)+(b*I*eps*(d1))*(PRUACC210)-(muND*PSDACCN*(ACC210ND/ND)*I*b*eps)-(muRD*PRDACC210)-(PRDACC210*(1/(8*365)))
+    #dACC210U.dt = del*(0.9)+(ACC210D/Tau)-(ACC210U*b*I*eps*d1)-(ACC210NU*muNU*I*b*eps) - (ACC210U*(1/(8*365)))
+    #dACC210D.dt = -(ACC210D/Tau)+ (ACC210U*b*I*eps*d1)-(ACC210ND*muND*I*b*eps) - (ACC210D*(1/(8*365)))
+    #dACC210NU.dt = del*(0.9)*(ACC210NU/ACC210U)+(ACC210ND/Tau)-(ACC210NU*b*I*eps*d1)-(ACC210NU*muNU*I*b*eps)- (ACC210NU*(1/(8*365)))
+    #dACC210ND.dt = -(ACC210ND/Tau)+ (ACC210NU*b*I*eps*d1)-(ACC210ND*muND*I*b*eps)- (ACC210ND*(1/(8*365)))
+    dPSUACCN2.dt = - (b*(PSU/N)*ws*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ND/Tau)*(PSDACCN2/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PSUACCN2/NU)#*(1-rho3)
+    dPSDACCN2.dt = - (b*(PSU/N)*ws*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) +(b*I*eps*NU*d1)*(PSUACCN2/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PSDACCN2/ND)- (muSD)*(PSDACCN2)#*(1-rho4)
+    dPRUACCN2.dt = - (b*(PRU/N)*h*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) +(ND/Tau)*(PRDACCN2/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PRUACCN2/NU)#*(1-rho3)
+    dPRDACCN2.dt = - (b*(PRU/N)*h*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (b*I*eps*NU*d1)*(PRUACCN2/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PRDACCN2/ND)#*(1-rho4)
+    dPSUACC.dt = - (b*(PSU/N)*ws*(SU*(meanstrainsS1+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(SU*(meanstrainsS1+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(b*I*eps*NU*rho1)*(PSUACCN2/NU)+(SD/Tau)*(PSDACC/SD)+(-(b*I*eps*SU*((d2*(meanstrainsS4+omega1))+(eta))) - (del*att*SU))*(PSUACC/SU)#*(1-rho3)
+    dPSDACC.dt = - (b*(PSU/N)*ws*(SD*(meanstrainsS2+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(SD*(meanstrainsS2+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ratr*b*I*eps*ND*rho2)*(PSDACCN2/ND)+(b*I*eps*SU*d2*(meanstrainsS4+omega1))*(PSUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*eta) - (del*att*SD))*(PSDACC/SD)- (muSD)*(PSDACC)#death rate included because of clearance by drugs rather than immune system
+    dPRUACC.dt = - (b*(PRU/N)*h*(SU*(meanstrainsS3+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SU*(meanstrainsS3+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (b*I*eps*NU*rho1)*(PRUACCN2/NU)+(SD/Tau)*(PRDACC/SD)+(-(b*I*eps*SU*((d2*(meanstrainsS4+omega1))+(eta))) - (del*att*SU))*(PRUACC/SU)#*(1-rho3)
+    dPRDACC.dt = - (b*(PRU/N)*h*(SD*(meanstrainsS4+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SD*(meanstrainsS4+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (ratr*b*I*eps*ND*rho2)*(PRDACCN2/ND)+(b*I*eps*SU*d2*(meanstrainsS4+omega1))*(PRUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*eta) - (del*att*SD))*(PRDACC/SD)#*(1-rho4)
+    dPSUACCA.dt = - (b*(PSU/N)*ws*(AU*(meanstrainsS5+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(AU*(meanstrainsS5+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(AD/Tau)*(PSDACCA/AD)+(b*I*eps*SU*eta)*(PSUACC/SU)-(b*I*eps*(meanstrainsS5)*AU*d3)*(PSUACCA/AU) - (del*att*AU)*(PSUACCA/AU)
+    dPSDACCA.dt = - (b*(PSU/N)*ws*(AD*(meanstrainsS6+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(AD*(meanstrainsS6+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ratr*b*I*eps*SD*eta)*(PSDACC/SD)+(b*I*eps*(omega1)*AU*d3)*(PSUACCA/AU)-(AD/Tau)*(PSDACCA/AD) - (del*att*AD)*(PSDACCA/AD) - (muSD)*(PSDACCA)
+    dPRUACCA.dt = - (b*(PRU/N)*h*(AU*(meanstrainsS5+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(AU*(meanstrainsS5+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(AD/Tau)*(PRDACCA/AD)+(b*I*eps*SU*eta)*(PRUACC/SU)-(b*I*eps*(omega1)*AU*d3)*(PRUACCA/AU)- (del*att*AU)*(PRUACCA/AU)
+    dPRDACCA.dt = - (b*(PRU/N)*h*(AD*(meanstrainsS6+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*h*(AD*(meanstrainsS6+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ratr*b*I*eps*SD*eta)*(PRDACC/SD)+(b*I*eps*(omega1)*AU*d3)*(PRUACCA/AU)-(AD/Tau)*(PRDACCA/AD) - (del*att*AD)*(PRDACCA/AD)
+    #dNUACC.dt = del+(ND/Tau)#-(b*I*eps*NU*(d1+rho1)) 
+    #dNDACC.dt = (b*I*eps*NU*d1)-(ND/Tau)#-(ratr*b*I*eps*ND*rho2)
+    return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,dbeta.dt,dPSUACCN2.dt,dPSDACCN2.dt,dPRUACCN2.dt,dPRDACCN2.dt,dPSUACC.dt,dPSDACC.dt,dPRUACC.dt,dPRDACC.dt,dPSUACCA.dt,dPSDACCA.dt,dPRUACCA.dt,dPRDACCA.dt)))
+  })
+}
+
+bavgnull = 1
+
+nstrains = 20
+
+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 = 1/100 #prob of appearance of new antigen
+omega2 = omega1 #prob of appearance of new antigen only in resistant strains
+
+
+L = 1
+
+ws = 1
+wr = 1-0.6
+#muSU = 1
+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 = 1
+d2 = d1#/6
+d3 = d1#/10
+#m = 3
+K=10
+p2 = list(bavgnull=bavgnull,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=40000,by=1) 
+
+
+N0 = c(500,500,500,500,333,333,333,333,333,333,bavgnull+bavgnull*L,200,200,200,200,200,200,200,200,200,200,200,200)#,100,100,70,70,150,150,150,150,500,500,500,500)#,500,500,500,500,150,150,150,150,200,200,200,200)
+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", ylim = c(0,max(out2[,5])),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")
+
+plot(out3[,1],out3[,12],xlim=c(0,2000))
+
+plot(out2[,1],out2[,17],type="l")
+plot(out2[,1],out2[,18],type="l")
+plot(out2[,1],out2[,19],type="l")
+plot(out2[,1],out2[,20],type="l")
+
+plot(out2[,1],out2[,21],type="l")
+plot(out2[,1],out2[,22],type="l")
+plot(out2[,1],out2[,23],type="l")
+plot(out2[,1],out2[,24],type="l")
+
+
+
+
+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 = matrix(data=0,nrow=6,ncol=366)
+out3eq2 = matrix(data=0,nrow=6,ncol=366)
+#out3eq3 = vector(mode = "numeric", length = 22)
+#out3eq4 = vector(mode = "numeric", length = 22)
+#out3eq5 = vector(mode="numeric", length=22)
+#out3eq6 = vector(mode="numeric", length=22)
+out3eq7 = matrix(data=0,nrow=6,ncol=366)
+m = matrix(data=0,nrow=6,ncol=366)
+Ivec = vector(mode="numeric",length = 6)
+mavg = vector(mode="numeric",length=6)
+favg = vector(mode="numeric",length=6)
+mnavg = vector(mode="numeric",length=19)
+#K = c(0,1,2,3,4,5,6,7,8,9,10)
+#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)
+bvalues = seq(from=-1.95,to=-1.45,by=0.025)
+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:6){
+  #del=delvalues[i]
+  bavgnull=10^(bvalues[i])
+  #wr=wrvalues[i]
+  #d1=dvalues[i]
+  #muSU = muSUvalues[i]*0.1
+  #rho3 = rho3values[i]
+  #omega1 = omegavalues[i]
+  wr=1-0.9
+  K = 10
+  p2 = list(bavgnull=bavgnull,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)
+  N0 = c(500,500,500,500,333,333,333,333,333,333,bavgnull+bavgnull*L,200,200,200,200,200,200,200,200,100,100,70,70)#,500,500,500,500,150,150,150,150,200,200,200,200)
+  out3 = ode(y=N0,times=t2,func=SEAI,parms=p2)
+  out3eq[i,] = out3[39635:40000,3]+out3[39635:40000,2] +out3[39635:40000,4] + out3[39635:40000,5]
+  #out3eq2[i,] = (out3[39635:40000,4]+out3[39635:40000,5])/(out3[39635:40000,2] + out3[39635:40000,3] + out3[39635:40000,4] + out3[39635:40000,5])
+  #out3eq3[i] = out3[399999,6]+out3[399999,7]
+  #out3eq4[i] = out3[399999,8]+out3[399999,9]
+  #out3eq5[i] = out3[399999,10]+out3[399999,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,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])
+  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,]))
+  #mnavg[i] = mean(1-exp(-(out3[39635:40000,17]+out3[39635:40000,18] +out3[39635:40000,19] + out3[39635:40000,20])/(out3[39635:40000,21]+out3[39635:40000,22])))
+  #ptwo = m[i]/(K-1)
+  #pthree = m[i]/(K-2)
+  #pfour = m[i]/(K-3)
+  #pfive = m[i]/(K-4)
+  #psix = m[i]/(K-5)
+  #pseven  = m[i]/(K-6)
+  #peight = m[i]/(K-7)
+  #pnine = m[i]/(K-8)
+  #pten = m[i]/(K-9)
+  #Ivec[i] = dnbinom(K,size=1,prob=pone) 
+  #favg[i] = mean(out3eq2[i,])
+}
+
+#plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant")
+#plot(bvalues,1-out4)#,ylab="Proportion Resistant")
+plot(bvalues,mavg)#,type="l")
+#plot(bvalues,mavg,type="l")
+#print(favg)
+print(mavg)
+plot(mavg,favg)
+
+#rho3values = c(1/20,1/50,1/100,1/250,1/400,1/500,1/750,1/1000,1/1250,1/1500)
+rho3values = seq(from = 0.005, to = 0.05, by = 0.005)
+mat6 = matrix(data=0,nrow=10,ncol=21)
+etavalues = seq(from=0.01,to=0.21,by=0.04)
+#mat77 = matrix(data=0,nrow=10,ncol=19)
+#mat7 = matrix(data=0,nrow=10,ncol=6)
+#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 5:6){
+  for (j in 1:21){
+    #d = dvalues[i]
+    #    A = Avalues[j]
+    #rho3=rho3values[i]
+    #rho4=rho3
+    rho3 = rho3values[i]
+    rho4=rho3
+    eta = 0.05
+    #wr = wrvalues[i]
+    wr=1-0.9
+    #muSU = muSUvalues[i]*0.01
+    bavgnull = 10^(bvalues[j])
+    #omega1 = omegavalues[j]
+    p2 = list(eta=eta,bavgnull=bavgnull,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)
+    N0 = c(500,500,500,500,333,333,333,333,333,333,bavgnull+bavgnull*L,200,200,200,200,200,200,200,200,200,200,200,200)#,100,100,70,70,150,150,150,150,500,500,500,500)#,500,500,500,500,150,150,150,150,200,200,200,200)
+    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] = mean((out5[39635:40000,4]+out5[39635:40000,5])/(out5[39635:40000,2] + out5[39635:40000,3] + out5[39635:40000,4] + out5[39635:40000,5]))
+    #mat7[i,j] = mean(1-exp(-(out5[39635:40000,2]+out5[39635:40000,3]+out5[39635:40000,4]+out5[39635:40000,5])/((out5[39635:40000,6]+out5[39635:40000,7]+out5[39635:40000,8]+out5[39635:40000,9]+out5[39635:40000,10]+out5[39635:40000,11]))))
+    #mat77[i,j] = mean(1-exp(-(out5[39635:40000,17]+out5[39635:40000,18] +out5[39635:40000,19] + out5[39635:40000,20])/(out5[39635:40000,21]+out5[39635:40000,22])))
+    #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(mat77)
+#print(mat8)
+#print(mat9)
+#print(mat10)
+#print(mat11)
+#print(mat12)
+
+plot(bvalues,mat6[1,],type="l",ylim=c(0,1))
+lines(bvalues,mat6[2,])
+lines(bvalues,mat6[3,])
+lines(bvalues,mat6[4,])
+lines(bvalues,mat6[5,])
+lines(bvalues,mat6[6,])
+lines(bvalues,mat6[7,])
+lines(bvalues,mat6[8,])
+lines(bvalues,mat6[9,])
+lines(bvalues,mat6[10,])
+
+plot(bvalues,mat77[9,],type="l",ylim=c(0,1))
+lines(bvalues,mat77[2,],col="RED")
+lines(bvalues,mat77[3,],col="GREEN")
+lines(bvalues,mat77[4,],col="BLUE")
+lines(bvalues,mat77[5,])
+lines(bvalues,mat77[6,])
+lines(bvalues,mat77[7,])
+lines(bvalues,mat77[8,])
+lines(bvalues,mat77[9,])
+lines(bvalues,mat77[10,])
+
+plot(mat77[1,],mat6[1,])
+plot(mat77[2,],mat6[2,])
+plot(mat77[3,],mat6[3,])
+plot(mat77[4,],mat6[4,])
+plot(mat77[5,],mat6[5,])
+plot(mat77[6,],mat6[6,])
+
+plot(mat77[9,2:7],mat6[9,])
+
+plot(mat77[1,9:12],mat6[1,9:12],type="l",ylim=c(0,1),xlim=c(0,1))
+lines(mat77[2,1:4],mat6[2,1:4])
+lines(mat77[3,5:8],mat6[3,5:8])
+
+plot(bvalues,mat7[1,],type="l",ylim=c(0,1))
+lines(bvalues,mat7[2,])
+lines(bvalues,mat7[3,])
+lines(bvalues,mat7[4,])
+lines(bvalues,mat7[5,])
+lines(bvalues,mat7[6,])
+lines(bvalues,mat7[7,])
+lines(bvalues,mat7[8,])
+lines(bvalues,mat7[9,])
+lines(bvalues,mat7[10,])
+
+plot(mat7[1,],mat6[1,],type="l",ylim=c(0,1))
+lines(mat7[2,],mat6[2,])
+lines(mat7[3,],mat6[3,])
+lines(mat7[4,],mat6[4,])
+lines(mat7[5,],mat6[5,])
+lines(mat7[6,],mat6[6,])
+lines(mat7[7,],mat6[7,])
+lines(mat7[8,],mat6[8,])
+lines(mat7[9,],mat6[9,])
+lines(mat7[10,],mat6[10,])
+
+
+
+
+
+
+
+mat7vec = c(mat7)
+mat7.df <- as.data.frame(c(t(mat7vec)))
+mat73.df <-cbind(mat7.df,vector50,vector62)
+
+mat77vec = c(mat77)
+mat77.df <-as.data.frame(c(t(mat77vec)))
+mat773.df <-cbind(mat77.df,vector50,vector62)
+#mat6[3,12] = 0
+#mat6[1,12] = 9.395332e-02
+#mat6[5,11] = 0
+#mat6[6,11] = 0
+mat6vec <- c(mat6)
+mat6.df <- as.data.frame(c(t(mat6vec)))
+vector50 <- c(rho3values,rho3values,rho3values,rho3values,rho3values,rho3values)#,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values)
+mat62.df <- cbind(mat6.df,vector50)
+
+vector60 = matrix(data=0,nrow=6,ncol=10)
+for(i in 1:6){
+  for(j in 1:10){
+    vector60[i,] <- rep(bvalues[i],10)
+  }
+}
+vector62 <- c(t(vector60))
+mat63.df <- cbind(mat62.df,vector62)
+
+rows <- rownames(mat63.df)
+
+columns <- colnames(mat63.df)
+log_transmission_intensity <- vector62
+nstrains_inverse <- vector50
+prevalence <- mat7vec
+frequency_resistant <- mat6vec
+naive_prevalence <- mat77vec
+ggplot(data=mat63.df,mapping=aes(x=log_transmission_intensity,y=frequency_resistant))+geom_line(mapping=aes(color=as.factor(nstrains_inverse)))
+
+ggplot(data=mat73.df,mapping=aes(x=log_transmission_intensity,y=prevalence))+geom_line(mapping=aes(color=as.factor(nstrains_inverse)))
+
+mixed.df <-cbind(mat63.df,mat7.df)
+
+ggplot(data=mixed.df,mapping=aes(x=prevalence,y=frequency_resistant))+geom_line(mapping=aes(color=as.factor(nstrains_inverse)))+geom_smooth()+geom_point()
+
+
+ggplot(data=mat773.df,mapping=aes(x=log_transmission_intensity,y=naive_prevalence))+geom_line(mapping=aes(color=as.factor(nstrains_inverse)))
+
+ggplot(data=mixed.df,mapping=aes(x=prevalence,y=frequency_resistant)) +geom_line(mapping=aes(color=as.factor(log_transmission_intensity)))+geom_smooth()
+
+mixed2.df <- cbind(mat63.df,mat77.df)
+
+ggplot(data=mixed.df,mapping=aes(x=nstrains_inverse,y=frequency_resistant)) +geom_line(mapping=aes(color=as.factor(log_transmission_intensity)))
+
+ggplot(data=mixed.df,mapping=aes(x=nstrains_inverse,y=frequency_resistant)) +geom_line(mapping=aes(color=as.factor(log_transmission_intensity)))
+
+ggplot(data=mixed.df,mapping=aes(x=nstrains_inverse,y=prevalence)) +geom_line(mapping=aes(color=as.factor(log_transmission_intensity)))
+
+ggplot(data=mixed2.df,mapping=aes(x=naive_prevalence,y=frequency_resistant)) + geom_point()+geom_smooth()
+
+vec600 = vector(mode="numeric",length=6)
+vec700 = vector(mode="numeric",length=6)
+vec770 = vector(mode="numeric",length=6)
+for (i in 1:6){
+  vec600[i] = mat6[i,i]
+  vec700[i] = mat7[i,i]
+  vec770[i] = mat77[i,i]
+}
+plot(vec700,vec600)
+
+require(ggplot2)
\ No newline at end of file