Skip to main content

Solutions: T2 – Building a model

🧠0% 🤖100%

These are solutions to the exercises at the end of T2 – Building a model. Many of the exercises are open-ended, so these are examples: there are often other good answers. This page runs on its own, from the setup below.

import starsim as ss

Exercise 1

A town of 2,000 people has had no measles for many years, and 70% of its people are immune. Five people arrive with measles. Using init and init_n, model the outbreak over a year. How many people are infected? Is the reproduction number at the start above or below 1?

With 30% susceptible, the reproduction number at the start is 15.2 × 0.3 ≈ 4.6 (printed as R = 4.53 at the start), well above 1, and almost all of the 600 susceptible people are infected. The run notes that with only 5 people infected at the start, whether an outbreak takes off is partly chance, which the ODE doesn’t show; running it many times with a stochastic solver (solver='gillespie', with ss.parallel(n_runs=...)) does.

S, I, R = ss.states('S I R')
infection = S >> I @ 1.9*I
recovery  = I >> R @ ss.days(8)
measles = ss.Disease(infection, recovery, init=dict(R=0.7), init_n=dict(I=5))

sim = ss.Sim(measles, total_pop=2_000, stop=ss.days(365))
sim.run()
print('Infections:', round(sim.results.measles.infection.sum()))
Sim: measles, 2,000 people, as an ODE, days 0 to 365
  measles: R0 = 15.2 (R = 4.53 at the start, with 30% susceptible), with 5 infectious at the start
  Note: "measles" starts with 5 infected, so whether it takes off is a matter of chance, which solver="ode" does not show; solver="gillespie" with ss.parallel(n_runs=...) does
  Done in 0.0051 s; measles: 588 infections; at most 270 infectious, on day 14
Infections: 588

Exercise 2

Using the final size relation, a model with R0 = 1.5 should infect about 58% of people. Choose β to give R0 = 1.5 with an infectious period of 8 days, run the model for long enough for the outbreak to end, and check.

β = R0 / the infectious period = 1.5 / 8 ≈ 0.19 per day. The outbreak is much slower than with R0 = 15, taking about two years to end.

S, I, R = ss.states('S I R')
infection = S >> I @ (1.5/8)*I
recovery  = I >> R @ ss.days(8)
measles = ss.Disease(infection, recovery, init_prev=0.001)

sim = ss.Sim(measles, total_pop=100_000, stop=ss.days(730))
sim.run()
print('Fraction infected:', round(sim.results.measles.R[-1] / 100_000, 2))
Sim: measles, 100,000 people, as an ODE, days 0 to 730
  measles: R0 = 1.5, with 100 infectious at the start
  Done in 0.0024 s; measles: 58,292 infections; at most 6,369 infectious, on day 91
Fraction infected: 0.58

Exercise 3

Change the model to density-dependent transmission (I.n), choosing β so that R0 is still 15 in a town of 100,000. Then run it in a town of 10,000 and of 1,000,000. How does the outbreak change with the size of the town, and how does this differ from the frequency-dependent model?

With density-dependent transmission, R0 = β N × the infectious period, so for R0 = 15.2 in a town of 100,000, β = 15.2 / (100,000 × 8) = 1.9 × 10⁻⁵ per infectious person per day. In a town of 10,000, R0 is then 1.5, and the slow outbreak infects about 60% of people; in a town of 1,000,000, R0 is 152, and almost everyone is infected within days. With frequency-dependent transmission, R0 and the fraction infected would be the same at every size. Each run prints its R0 with verbose=1, and sim.diseases.measles.R0() gives it.

beta = 1.9e-5 # Per infectious person per day
for pop in [10_000, 100_000, 1_000_000]:
    S, I, R = ss.states('S I R')
    infection = S >> I @ beta*I.n
    recovery  = I >> R @ ss.days(8)
    measles = ss.Disease(infection, recovery, init_n=dict(I=10))
    sim = ss.Sim(measles, total_pop=pop, stop=ss.days(365), verbose=0)
    sim.run()
    res = sim.results.measles
    print('Town of', pop, ' fraction infected:', round(res.R[-1] / pop, 2), ' peak infectious:', round(res.I.max()))
Town of 10000  fraction infected: 0.6  peak infectious: 673
Town of 100000  fraction infected: 1.0  peak infectious: 75229
Town of 1000000  fraction infected: 1.0  peak infectious: 951749