using Gadfly
using DataFrames
using Setfield
set_default_plot_size(18cm, 12cm)
From https://en.wikipedia.org/wiki/Mathematical_modelling_of_infectious_disease#The_SIR_model
${\displaystyle S(t)}$ is used to represent the number of individuals not yet infected with the disease at time t, or those susceptible to the disease.
${\displaystyle I(t)}$ denotes the number of individuals who have been infected with the disease and are capable of spreading the disease to those in the susceptible category.
${\displaystyle R(t)}$ is the compartment used for those individuals who have been infected and then removed from the disease, either due to immunization or due to death. Those in this category are not able to be infected again or to transmit the infection to others.
${\displaystyle \beta }$ is the rate of contracting the disease
$\gamma$ represents the mean recovery/death rate, or ${\displaystyle {1/{\gamma}}}$ the mean infective period
struct CVconf
N
beta
gamma
end
struct CVstate
t
S
R
I
end
CVstate(conf, I) = CVstate(0, conf.N, 0, I)
function evolve_step(old::CVstate, conf)
dS = -conf.beta * old.S * old.I / conf.N
dR = conf.gamma * old.I
dI = - dS - dR
CVstate(old.t + 1, old.S + dS, old.R + dR, old.I + dI)
end
evolve_step (generic function with 1 method)
function evolve(state, conf, time)
r = []
for i in 1:time
state = evolve_step(state, conf)
push!(r, state)
end
r
end
evolve (generic function with 1 method)
function evolve_seq(state, conf, seq)
r = []
for s in seq
time, beta = s
sconf = @set conf.beta = beta
append!(r, evolve(sconf, state, time))
state = r[end]
end
r
end
evolve_seq (generic function with 1 method)
function plot_evolution(data, args...)
plot(stack(DataFrame(data)), x=:t, y=:value, color=:variable, Geom.line,
Scale.y_log10,
Guide.xticks(ticks=0:7:data[end].t), Guide.yticks(ticks=0:8),
Coord.cartesian(ymin=2),
Guide.ylabel("people"), Guide.xlabel("time / days"),
args...)
end
plot_evolution (generic function with 1 method)
conf = CVconf(8e6, 0.3, 1/7)
r = evolve(CVstate(conf, 100), conf, 284);
#plot_evolution(r)
struct CVDetailConf
N
beta
k_severe
gamma
k_recover
k_dead
end
struct CVDetailState
t
S
Rdead
Ralive
Imild
Isevere
end
CVDetailState(conf, I) = CVDetailState(0, conf.N, 0, 0, I, 0)
function evolve_step(old::CVDetailState, conf)
# S -beta*[Imild]-> Imild
# Imild -k_severe-> Isevere
# Imild -gamma-> Ralive
# Isevere -gamma*k_dead-> Rdead
# Isevere -gamma*(1-k_dead)-> Ralive
dS = -conf.beta * old.S * old.Imild / conf.N
dRmild_alive = conf.gamma * old.Imild
dRsevere_alive = conf.k_recover * old.Isevere
dRsevere_dead = conf.k_dead * old.Isevere
dImild_severe = conf.k_severe * old.Imild
dIsevere = dImild_severe - dRsevere_alive - dRsevere_dead
dImild = - dS - dRmild_alive - dImild_severe
CVDetailState(old.t + 1,
old.S + dS,
old.Rdead + dRsevere_dead,
old.Ralive + dRmild_alive + dRsevere_alive,
old.Imild + dImild,
old.Isevere + dIsevere
)
end
evolve_step (generic function with 2 methods)
The model parameters were chosen to match the few numbers we know fairly certain, from Italy:
The following plots try to estimate the number of total asymptomatic or low symptomatic patients.
The first plot chooses a low severity rate, which leads to a high amount of undiagnosed people, but also leads to a low case fatality rate.
The second plot chooses a high severity rate, which leads to relatively few undiagnosed people; however, the result is a high case fatality rate.
Both models are wrong, but how wrong, and why? Please let me know if you have ideas.
conf_low = CVDetailConf(60e6, 0.6, 0.005, 1/7, 1/1.5, 1/3)
r = evolve(CVDetailState(conf_low, 10), conf_low, 200)
plot_evolution(r, Guide.title("low severity rate, high # of undiagnosed cases"))
┌ Warning: `getindex(df::DataFrame, col_ind::ColumnIndex)` is deprecated, use `df[!, col_ind]` instead. │ caller = evalmapping(::DataFrame, ::Symbol) at dataframes.jl:96 └ @ Gadfly /home/corecode/.julia/packages/Gadfly/1wgcD/src/dataframes.jl:96
conf_high = CVDetailConf(60e6, 0.6, 0.15, 1/7, 1/1.5, 1/4)
r = evolve(CVDetailState(conf_high, 10), conf_high, 200)
plot_evolution(r, Guide.title("high severity rate"))