In this activity, we are trying to simulate an outbreak of a new infectious disease that our population of 10^6people has not been exposed to before. This means that we are starting with a single case, everyone else is susceptible to the disease, and no one is yet immune or recovered. This can for example reflect a situation where an infected person introduces a new disease into a geographically isolated population, like on an island, or even when an infections "spill over" from other animals into a human population. In terms of the initial conditions for our model, we can define: "S=10^(6) -1=999999," I = 0and R = 0. NULL

Remember, the differential equations for the simple SIR model look like this:

dS/dt = `λS`*dI/dt and `λS`*dI/dt = `λS`-I*gamma, dR/dt = I*gamma

Initial number of people in each compartment
S = 10^6-1",I=0  "and R = 0.

NULL

Parameters:

gamma = .1*recovery*rate*beta and .1*recovery*rate*beta = .4*the*daily*infection*rate

restart; with(plots); _local(gamma)

sys := diff(s(t), t) = -lambda*s(t), diff(i(t), t) = lambda*s(t)-gamma*i(t), diff(r(t), t) = gamma*i(t)

diff(s(t), t) = -lambda*s(t), diff(i(t), t) = lambda*s(t)-gamma*i(t), diff(r(t), t) = gamma*i(t)

(1)

ic := s(0) = s__0, i(0) = i__0, r(0) = r__0

gamma := .1; beta := .4; n := 10^6

.1

 

.4

 

1000000

(2)

lambda := beta*i(t)/n

s__0, i__0, r__0 := 10^6-1, 1, 0

NULL

sols := dsolve({ic, sys}, numeric, output = listprocedure)

display([odeplot(sols, [t, s(t)], 0 .. 100, color = red), odeplot(sols, [t, i(t)], 0 .. 100, color = blue), odeplot(sols, [t, r(t)], 0 .. 100, color = green)], labels = ["Time [day]", "Population"], labeldirections = [horizontal, vertical], legend = ["Susceptible", "Infected", "Recovered"], legendstyle = [location = right])

 

Remember that in a simple homogenous SIR model, `R__eff  `is directly related to the proportion of the population that is susceptible:

R__eff = R__0*S/N

Reff := proc (t) options operator, arrow; beta*s(t)/(gamma*n) end proc

odeplot(sols, [[t, Reff(t)]], t = 0 .. 100, size = [500, 300], labels = ["Time [day]", "Reff"], labeldirections = [horizontal, vertical])

 

The effective reproduction number is highest when everyone is susceptible: at the beginning, `R__eff  ` = R__0. At this point in our example, every infected cases causes an average of 4 secondary infections. Over the course of the epidemic, `R__eff  ` declines in proportion to susceptibility.

The peak of the epidemic happens when `R__eff  ` goes down to 1 (in the example here, after 50 days). As `R__eff  `decreases further below 1, the epidemic prevalence goes into decline. This is exactly what you would expect, given your understanding of the meaning of `R__eff  ` once the epidemic reaches the point where every infected case cannot cause at least one more infected case (that is, when `R__eff  ` < 1), the epidemic cannot sustain itself and comes to an end.

susceptible := eval(s(t), sols); infected := eval(i(t), sols); recovered := eval(r(t), sols)

susceptible(51)

HFloat(219673.04834159758)

(3)

infected(51)

HFloat(401423.4112878752)

(4)

recovered(51)

HFloat(378903.54037052736)

(5)

Reffe := proc (t) options operator, arrow; beta*susceptible(t)/(gamma*n) end proc

proc (t) options operator, arrow; beta*susceptible(t)/(gamma*n) end proc

(6)

Reffe(51)

HFloat(0.8786921933663903)

(7)

Prevalence is simply the value of Iat a given point in time. Now we can see that the incidence is the number of new cases arriving in the I compartment in a given interval of time. The way we represent this mathematically is by taking the integral of new cases over a given duration.

For example, if we wanted to calculate the incidence from day 7 to 14,

int(`&lambda;S`(t), t = 7 .. 14)

lamda := proc (t) options operator, arrow; beta*infected(t)/n end proc

proc (t) options operator, arrow; beta*infected(t)/n end proc

(8)

inflow := proc (t) options operator, arrow; lamda(t)*susceptible(t) end proc

proc (t) options operator, arrow; lamda(t)*susceptible(t) end proc

(9)

int(inflow(t), t = 7 .. 14)

HFloat(78.01804723222038)

(10)

incidence_plot := plot(inflow(t), t = 0 .. 14, color = orange, labels = ["Time (days)", "Incidence Rate"], labeldirections = [horizontal, vertical], title = "Incidence Rate between t=7 and t=14")

 

s, i, r := eval(s(t), sols), eval(i(t), sols), eval(r(t), sols); T := 100; dataArr := Array(-1 .. T, 1 .. 4); dataArr[-1, () .. ()] := `<,>`("Day", "Susceptible", "Infected", "Recovered")


Assign all the subsequent rows

for t from 0 to T do dataArr[t, () .. ()] := `~`[round](`<,>`(t, s(t), i(t), r(t))) end do

 

Tabulate through the DocumentTools

DocumentTools:-Tabulate(dataArr, alignment = left, width = 50, fillcolor = (proc (A, n, m) options operator, arrow; ifelse(n = 1, "DeepSkyBlue", "LightBlue") end proc))

Download dynamics_of_novel_disease_outbreak.mw


Please Wait...