Skip to main content

Hospitalization and death

🧠0% 🤖100%

Two different mechanisms decide what happens to people after a state, and a model needs the right one. In a split on entry, the outcome is decided when people enter the state: here, 10% of infectious people will go to hospital and the rest will recover, each after their infectious period. With competing hazards, two transitions out of a state run at the same time and whichever happens first wins: here, people in hospital are discharged after 10 days on average, but die at a rate of 0.02 per day while they are there, so the fraction who die depends on both rates. This page puts both in one model, a respiratory infection with hospitalization, and checks the fraction who die. Then it makes the chance of hospitalization depend on age, written as a Python function, and runs the model both as compartments and as agents.

The model

Parameter Value Meaning
β 0.3 per day Rate of infection
Infectious period 8 days Mean time in I (exponential), before recovery or hospitalization
Hospitalized 10% Fraction of infectious people who go to hospital, decided on entry to I
Stay in hospital 10 days Mean time to discharge
Death in hospital 0.02 per day Rate of death while in hospital, competing with discharge
R0 2.4 β × the infectious period
Population 100,000
Infectious at the start 0.1%

The model is written as a function of the chance of hospitalization, so that the same model can be made again with a chance that depends on age, below:

import starsim as ss

def make_covid(hospitalized=0.1):
    """ An SIR with hospitalization and death; hospitalized is the chance of going to hospital """
    S, I, R, H, D = ss.states('S I R H D')
    infection = S >> I @ 0.3*I
    course    = ss.Transition(I, {R: None, H: hospitalized}, dur=ss.days(8))  # Split on entry: 10% to H, the rest to R
    discharge = H >> R @ ss.days(10)                                            # Competing hazards: discharge ...
    death     = H >> D @ ss.perday(0.02)                                        # ... and death
    return ss.Disease(infection, course, discharge, death, init=dict(I=0.001), name='covid')

covid = make_covid()

In {R: None, H: 0.1}, None stands for the rest (here, 90%). The diagram, from diagram():

covid.diagram()

And the equations, in which the split appears as the fractions (1 - 0.1) and 0.1 of the flow out of I, and the competing hazards as two flows out of H:

covid.equations()
dS/dt = -beta * S * I / N
dI/dt =  beta * S * I / N - I / dur_I
dR/dt =  (1 - 0.1) * I / dur_I + H / dur_H
dH/dt =  0.1 * I / dur_I - rate_death * H - H / dur_H
dD/dt =  rate_death * H

A rate given as a number is named by its transition, so the rate of death appears as a parameter called rate_death, and a duration by the state it is the time in (dur_I and dur_H). To give it a name of your own, write it as ss.par(ss.perday(0.02), name='mu') (see Parameters).

Result

sim = ss.Sim(covid, total_pop=100_000, stop=ss.days(200))
sim.run()
sim.plot('covid')
Sim: covid, 100,000 people, as an ODE, days 0 to 200
  covid: R0 = 2.4, with 100 infectious at the start
  Done in 0.0025 s; covid: 87,777 infections; at most 21,892 infectious, on day 41

Of the people admitted to hospital, the fraction who die is the rate of death over the sum of the two rates out of H: 0.02 / (0.02 + 1/10) = 1/6. The flow course.H is the number admitted to hospital:

admitted = sim.results.covid['course.H'].sum()
deaths = sim.results.covid.death.sum()
print('Admitted:', round(admitted), 'Died:', round(deaths), 'Fraction:', round(deaths/admitted, 3))
Admitted: 8788 Died: 1465 Fraction: 0.167

Advanced: severity by age, as a function

Some parts of a model are easiest to write as a function. Here the chance that an infection is severe enough for hospital depends on age: 10% under 50, and 40% from 50 up. Rather than stratify the disease by age and give a value for each group, the chance is a Python function of the sim, which reads people’s ages. Starsim tries the function on the model’s age groups, and since it works on them, the model still runs as compartments; run as agents, the same function is evaluated for each person.

def severity(sim):
    """ The chance of hospitalization: 10%, plus 30% from age 50 """
    return 0.1 + 0.3*(sim.people.age >= 50)

covid_age = make_covid(hospitalized=severity)

The sim is split into two age groups at 50 (stratify=dict(age=[0, 50])), so that the function has ages to work on in a compartmental run. The same model is run as an ODE and as agents, in a population of 10,000:

sims = {
    'ODE': ss.Sim(covid_age, total_pop=10_000, stop=ss.days(200), stratify=dict(age=[0, 50]), verbose=0),
    'Agents': ss.Sim(covid_age, total_pop=10_000, stop=ss.days(200), solver='abm', verbose=0),
}
msim = ss.parallel(sims)
msim.plot('covid')

The two agree, apart from chance in the agent run. Of about 8,800 infections, about 1,500 people are admitted to hospital (covid_course.H): about one in six, between the 10% of those under 50 and the 40% of those over, since about a quarter of people are 50 or over (the population’s ages are guessed as the world’s in 2024, as describe() says). Of those admitted, one in six die, as before, since the competing hazards out of H haven’t changed:

msim.summarize().loc[['covid_infection', 'covid_course.H', 'covid_death']]
value
result scenario
covid_infection ODE 8,778
Agents 8,795
covid_course.H ODE 1,547
Agents 1,529
covid_death ODE 258
Agents 247

A function that needs each person (one that takes uids, the people it applies to, as well as the sim) can’t be run as compartments, so Starsim runs a model with one as agents, and describe() says why.

Going further