sir-with-variant-heterogeneous.sage

## Consider a population infected by a virus which has two forms, an
## ancestral one and a variant (with perfect cross-immunity) having a
## different basic reproduction number and possibly a different
## recovery rate.  We simulate the evolution of the epidemic possibly
## with heterogeneous susceptibility and/or epidemic quenching (viz.,
## mitigation measures which depend on number of infected).

## s  = susceptible
## i1 = infected by ancestral form
## i2 = infected by variant
## r  = recovered (immune to both)
## sx2= not having been specifically infeced by variant (note: included in s)

var('s i1 i2 r sx2 t')

# Reproduction number of the ancestral form
rep1 = 1.0
# Reproduction number of the variant for the subpopulation
rep2 = 1.4

# Recovery rate for variant
gamma2 = 1.0

# Initial proportion infected
eps = 0.0075
# Proportion OF infected initially by variant
vprop = 1e-2

# Quenching parameter (use 0 for standard SIR)
quench = 0.0

# Heterogeneity coefficient (1 = standard SIR, 2 = Gamma-distributed)
hetcoef = 1.0
# Extra heterogeneity for variant
hetcoefx2 = 0.0

beta1 = rep1
beta2 = rep2*gamma2
sdot = -(beta1*i1+beta2*i2*sx2^hetcoefx2)*s^hetcoef*(1-quench*(i1+i2))
i1dot = beta1*s^hetcoef*i1*(1-quench*(i1+i2))-i1
infect2 = beta2*s^hetcoef*i2*sx2^hetcoefx2*(1-quench*(i1+i2))
i2dot = infect2-gamma2*i2
rdot = i1+gamma2*i2
p = desolve_system_rk4([sdot, i1dot, i2dot, rdot, -infect2], [s,i1,i2,r,sx2], ics=[0, 1-eps, eps*(1-vprop), eps*vprop, 0, 1-eps*(1-vprop)], ivar=t, end_points=[0, 50], step=0.1)

g1 = line([(tmp[0], tmp[1]) for tmp in p], color="green") + line([(tmp[0], tmp[2]+tmp[3]) for tmp in p], color="red") + line([(tmp[0], tmp[4]) for tmp in p], color="blue") + line([(tmp[0], tmp[3]/(tmp[2]+tmp[3])) for tmp in p], color="purple", linestyle="dashed")
g2 = line([(tmp[0], tmp[2]+tmp[3]) for tmp in p], color="red", ymin=0, ymax=0.1) + line([(tmp[0], tmp[2]) for tmp in p], color="orange", ymin=0, ymax=0.1) + line([(tmp[0], tmp[3]) for tmp in p], color="purple", ymin=0, ymax=0.025)
g3 = line([(tmp[0], tmp[1]) for tmp in p], color="green", ymin=1e-4, ymax=1, scale="semilogy") + line([(tmp[0], tmp[2]+tmp[3]) for tmp in p], color="red", ymin=1e-4, ymax=1, scale="semilogy") + line([(tmp[0], tmp[4]) for tmp in p], color="blue", ymin=1e-4, ymax=1, scale="semilogy") + line([(tmp[0], tmp[3]/(tmp[2]+tmp[3])) for tmp in p], color="purple", linestyle="dashed", ymin=1e-4, ymax=1, scale="semilogy")
g4 = line([(tmp[0], 1) for tmp in p], color="gray", linestyle="dashdot", ymin=0.5, ymax=1.5) + line([(tmp[0], (1+i1dot/i1).subs({s:tmp[1], i1:tmp[2], i2:tmp[3], sx2:tmp[5]})) for tmp in p], color=Color("orange").darker(), thickness=0.5, ymin=0.5, ymax=1.5) + line([(tmp[0], (1+i2dot/(gamma2*i2)).subs({s:tmp[1], i1:tmp[2], i2:tmp[3], sx2:tmp[5]})) for tmp in p], color=Color("purple").darker(), thickness=0.5, ymin=0.5, ymax=1.5) + line([(tmp[0], (1+(i1dot+i2dot)/(i1+i2)).subs({s:tmp[1], i1:tmp[2], i2:tmp[3], sx2:tmp[5]})) for tmp in p], color="black", linestyle="dotted", ymin=0.5, ymax=1.5) + line([(tmp[0], (1+(i1dot+i2dot)/rdot).subs({s:tmp[1], i1:tmp[2], i2:tmp[3], sx2:tmp[5]})) for tmp in p], color="black", ymin=0.5, ymax=1.5)
graphics_array([[g1, g2], [g3, g4]]).show(dpi=192)