sir-with-variant.sage

## Consider a population infected by a virus which has two forms, an
## ancestral one and a variant (with perfect cross-immunity): the
## population is equally susceptible to the ancestral form, but a
## certain subpopulation is much more highly susceptible to the
## variant.

## s1 = "normal" susceptible
## s2 = more susceptible to variant
## i1 = infected by ancestral form
## i2 = infected by variant
## r  = recovered (immune to both)

var('s1 s2 i1 i2 r t')

# Proportion more susceptible to the variant:
prop2 = 1.0
# Reproduction number of the ancestral form (in a wholly susceptible population!)
rep1 = 1.5
# Reproduction number of the variant for the subpopulation
rep2 = 2.0

# Initial proportion infected
eps = 0.001
# Proportion OF infected initially by variant
vprop = 0.001

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

s1dot = -(rep1*i1+rep1*i2)*s1*(1-quench*(i1+i2))
s2dot = -(rep1*i1+rep2*i2)*s2*(1-quench*(i1+i2))
i1dot = (rep1*s1+rep1*s2)*i1*(1-quench*(i1+i2))-i1
i2dot = (rep1*s1+rep2*s2)*i2*(1-quench*(i1+i2))-i2
rdot = i1+i2
p = desolve_system_rk4([s1dot, s2dot, i1dot, i2dot, rdot], [s1,s2,i1,i2,r], ics=[0, (1-prop2)*(1-eps), prop2*(1-eps), eps*(1-vprop), eps*vprop, 0], ivar=t, end_points=[0, 30])

g1 = line([(tmp[0], tmp[1]+tmp[2]) for tmp in p], color="green") + line([(tmp[0], tmp[3]+tmp[4]) for tmp in p], color="red") + line([(tmp[0], tmp[5]) for tmp in p], color="blue") + line([(tmp[0], tmp[4]/(tmp[3]+tmp[4])) for tmp in p], color="purple", linestyle="dashed")
g2 = line([(tmp[0], tmp[3]+tmp[4]) for tmp in p], color="red", ymin=0, ymax=0.1) + line([(tmp[0], tmp[3]) for tmp in p], color="orange", ymin=0, ymax=0.1) + line([(tmp[0], tmp[4]) for tmp in p], color="purple", ymin=0, ymax=0.1)
g3 = line([(tmp[0], tmp[1]+tmp[2]) for tmp in p], color="green", ymin=1e-4, ymax=1, scale="semilogy") + line([(tmp[0], tmp[3]+tmp[4]) for tmp in p], color="red", ymin=1e-4, ymax=1, scale="semilogy") + line([(tmp[0], tmp[5]) for tmp in p], color="blue", ymin=1e-4, ymax=1, scale="semilogy") + line([(tmp[0], tmp[4]/(tmp[3]+tmp[4])) for tmp in p], color="purple", linestyle="dashed", ymin=1e-4, ymax=1, scale="semilogy")
g4 = line([(tmp[0], (1+(i1dot+i2dot)/(i1+i2)).subs({s1:tmp[1], s2:tmp[2], i1:tmp[3], i2:tmp[4]})) for tmp in p], color="black", ymin=0.0, ymax=2.0)
graphics_array([[g1, g2], [g3, g4]]).show(dpi=192)