Last active
March 14, 2021 17:52
-
-
Save Gro-Tsen/4eacb2217f0fda0042a9eeab2a503e22 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## 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) | |
var('s i1 i2 r 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 | |
beta1 = rep1 | |
beta2 = rep2*gamma2 | |
sdot = -(beta1*i1+beta2*i2)*s^hetcoef*(1-quench*(i1+i2)) | |
i1dot = beta1*s^hetcoef*i1*(1-quench*(i1+i2))-i1 | |
i2dot = beta2*s^hetcoef*i2*(1-quench*(i1+i2))-gamma2*i2 | |
rdot = i1+gamma2*i2 | |
p = desolve_system_rk4([sdot, i1dot, i2dot, rdot], [s,i1,i2,r], ics=[0, 1-eps, eps*(1-vprop), eps*vprop, 0], 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]})) 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]})) 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]})) 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]})) for tmp in p], color="black", ymin=0.5, ymax=1.5) | |
graphics_array([[g1, g2], [g3, g4]]).show(dpi=192) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment