COVID-19 modeling: multiple populations: model description and Wuhan calibration

the effect of poorly-behaved neighbors

13 min readApr 10, 2020

introduction

Here I describe my model for treating interacting groups with Alison Hill’s SEIR model for COVID-19. It’s a bit heavy on details, so I’ve segregated the model description here, and will present some results in a following article. I think of this as a “red state, blue state” model, but it’s generally applicable to any situation where there’s inhomogeneity in people’s response to an outbreak.

This is likely not original work: I’ve done no literature search, and I’m not an epidemiologist. I did this as an exercise for myself. I’m sure all of this stuff was worked out a hundred years ago.

basic model

Alison Hill published a COVID-19 model which I adopted to PERL. The model, since enhanced, is an SEIR-type model, with the following stages, one of which I added:

S: susceptible: Not yet infected, with no immunity
E: exposed: Infected, showing no symptoms, and not yet infectious. This period was assumed to last 3 day, until patient proceeds to P.
P: pre-symptomatic: Still no symptoms, but now infectious. I assumed based on published data this period lasts 2.5 days, at which point the patient either recovers or proceeds to I1.
I1: mild symptomatic: Symptoms have begun, but hospitalization may not be required yet except possibly for quarantine. This period lasts 4 days, on average, according to the original parameters, after which the patient proceeds to either I2 or recovers. In my fitting, I assumed infectiousness was half what it was during the P period, to remain consistent with reported doubling times and R₀ numbers.
I2: serious symptomatic: No longer infectious (due to no longer interacting with others, and due to less active virus). This period lasts 4 days, after which patient either recovers or moves to I3.
I3: critical symptomatic: The patient is gravely ill, and not infectious since isolated in the hospital and viral levels may be much lower. This period lasts on average 7.7 days, after which the patient either dies (31%) or recovers (69%)

The differential equations are the following, where I added my P phase, as well keeping count of A, the number of cases which are infected and infectious but recover without symptoms:

∂S/∂t = –(β₀ P + β₁ I1+ β₂ I2+ β₃ I3) S
∂E/∂t = (β₀ P + β₁ I1+ β₂ I2+ β₃ I3) S — p₀ E
∂P/∂t = p₀ E — (σ + γ₀) P
∂I1/∂t = σ P— (γ₁ + p₁) I1
∂I2/∂t = p₁ I1— (γ₂ + p₂) I2
∂I3/∂t = p₂ I2 — (γ₃ + μ) I3
∂R/∂t = γ₀ P + γ₁ I1+ γ₂ I2 + γ₃ I3
∂A/∂t = γ₀ P
∂D/∂t = μ I3

As is clear from the equations of these segments except for S has a duration which is exponentially distributed, not constant, so while the listed durations are averages, the actual durations in a specific case can be shorter or longer. I tuned the coefficients to give what seemed like “reasonable” means based on what I read, but don’t take them as definitive. I based these estimates on various references, but I compare them to a very recently published CDC paper in the next section, and they seem reasonable.

The boldest thing I did was to assign double the infectiousness to the P period as I did to the I1 period. I did this to keep R₀ in the range 2.5 to 3, which has been widely reported, while allowing a doubling time of 3 days, yet with a mean time to infectiousness from infection of 2.5 days. With such a long lag relative to the doubling time, either patients are infecting a lot more than 2 others, or they infect people over a relatively short time (the serial interval).

The CDC paper (Sanche, et al.)

This recent study at CDC used a much higher value of R₀ = 5.7, in part because of their relatively long serial interval. If the mean time to infection is longer, each person needs to infect more others, to get the same doubling time. For example, if people infect people after one doubling time, they must be infecting 2 others before recovering. But if they’re infecting them after two doubling times, they must be infecting 4 others. Etc. This would have substantial implications on the magnitude of interventions needed to get R₀ < 1 and therefore to attain a reduction in cases over time.

This CDC study also has some useful time interval numbers:

infection to symptoms: 4.5 days (I used 5.5 days)
symptoms to hospitalization (serious infection): 5.5 days (I used 4.0 days)
hospitalization to death: 11 days (I used 11.7 days)
infection to death: 21 days (I used 21 days)

So my parameters appear to be in the ballpark of the ones used by CDC, with the exception I used a higher β in the pre-infectious period to reduce the effective serial interval and allow for a lower R₀

The main variable in this model are the β values, which represent the likelihood per unit time an infectious person infects a susceptible person. This depends on three factors: how much susceptible people are in public, how much infectious people are in public, and how likely it is that an infectious person will infect susceptible people in proximity.

interventions

the β values are the most obvious parameter to change during interventions. Those who are susceptible or pre-symptomatic are less likely to encounter others, and those who have symptoms and are infectious are more likely to self-quarantine at the appearance of symptoms. So to implement sanctions, I reduce the β values, keeping my ratio of 2:1 for the β’s associated with the P : I1 phases.

Wuhan data

I adjusted parameters of the model to fit data from Wuhan, which has been brought into question. But they were available, and allowed for experimentation with the model. Here are the results:

There’s three phases here, so I set my code to allow for various phases of β multiplication factors:

100% before interventions,
34% after moderate interventions,
8% with strong interventions.

Here are the baseline parameters I used in this model:

$par{mu}      = 0.04;
$par{beta0}   = 0.854;
$par{beta1}   = 0.427;
$par{beta2}   = 0.00;
$par{beta3}   = 0.00;
$par{gamma0}  = 0.2;
$par{gamma1}  = 0.200;
$par{gamma2}  = 0.188;
$par{gamma3}  = 0.090;
$par{p0}      = 0.333;
$par{sigma}   = 0.4;
$par{p1}      = 0.050;
$par{p2}      = 0.062;

The model results are shown with solid lines, with the bands representing +/- 2σ variation from Poisson statistics, assuming each positive result is independent (no clustered tests). The model captures the reported data well (maybe too well).

Some numbers from this simulation:

R0                    = 2.5620
preinfectious time    = 3.0030
presymptomatic time   = 2.5000
incubation time       = 5.5030
asymptomatic fraction = 0.3333
mild fraction         = 0.8000
mild duration         = 4.0000
severe fraction       = 0.2000
severe duration       = 4.0000
critical fraction     = 0.0496
critical duration     = 7.6923
symptoms to death     = 15.6923
infections to death   = 21.1953
case fatality         = 0.0102
symptomatic transmission fraction = 0.3902group 0 doubling time = 3.6425
group 0: exponential regime:
E   = 53.878%
P   = 22.702%
I1  = 20.624%
I2  = 2.342%
I3  = 0.453%

One important factor here: although the exposed and pre-symptomatic infectious phases are relatively short, during exponential growth (solved with an eigenvalue problem), they are responsible for 76.6% of the active cases, even with the modest doubling time of 3.6 days in these data. This is because exponential growth skews the case distribution earlier in the cycle.

multiple groups

An over-simplification of this model is it assumes a population is homegeneous, meaning everyone behaves the same. This isn’t realistic. In real life, some people are sheltered, working from home, venturing forth only for food shopping (if that). Others are flaunting restrictions and crowding the beaches. Others are essential workers who are doing their best in grocery stores. Others provide emergency services and are attending to sick, infectious patients.

But the real motivation for me to add this is the red state — blue state problem. Governors of blue states, for example Gavin Newsom of California, tended to be relatively quick to demand social isolation, while red state governors, like DeSantis of Florida, a buffoon of the highest order, were very late.

So I assume multiple groups, each with their own β values. In particular I assume the same starting β values, but then different group adapt interventions with potentially different strength and at different times.

group interactions

With non-interacting groups, the result is trivial: just model each group separately. This would be applicable for different continents where there is no longer travel between the continents, no virus transmission through shipments, etc.

The Chris Murray Model (University of Washington) treats each state as a separate, non-interacting group. It assumes if one state enacts strong interventions, it will have a peak, while if another state delays enacting interventions, that state will have a delayed peak. It however assumes that members of the latter state will not infect anyone from the former. This is a simplification and can miss the importance of enacting national, rather than local, measures.

Here I am trying to allow for the case where there are different groups which may mostly interact among themselves, but do have some interaction with each other, to examine the extent to which populations with higher β values can retard the recovery of regions with lower β values.

β matrix for intergroup infections

With Alison Hill’s model, β is a one-dimensional matrix, sort of a vector, describing how rapidly each potentially infectious phase causes susceptible people to become exposed. In the model, with my addition of a β₀ for the P-phase, there are four β values: β₀, β₁, β₂, and β₃. Following Alison’s defaults, I set β₂ =β₃ = 0, but the model supports positive values for these. Large values for these would increase R₀, however.

With multiple groups, however, there are separate β values for each phase infecting susceptible members of each of the other groups. So if I enumerate groups as a, b, c … then β₀ab might be the rate at which members of the P population of group a infect the S population of group b. So if I have N groups, this results in 4N² coefficients, which is a lot of β values.

β matrix: interventions

I want to simplify the number of coefficients, based on some reasonable assumptions.

First, I assume each group has the same basic coefficients for themselves, and the differences between groups are determined by a scale factor for that group at that particular time. So each group starts with the same values of β₀, β₁, β₂, and β₃. When a group practices an interventions, its β values may scale by a factor. Other groups may scale their β values different amounts at different times.

But I still need to interaction components. If I have two groups with different β values, how do I “guess” what the interaction components are?

β matrix: mixing factor

I start with the assertion that the strength of interaction between groups is characterized by a scale factor, mab. So if mab = 1, then groups a and b maximally interact as if they were subgroups of the same group, while if mab = 0, then the groups fully isolate from each other (but may still be connected via a third group c, where mac > 0 and mcb > 0, for example).

To do this, I have two trivial limits. These aren’t assumption, but are necessary conditions:

if a value mab = 0, then the β matrix elements are the same as each group on its own, for within-group infections, and zero for inter-group infections.
if a value mab = 1, and each of the two groups has the same β values, and the net result is the same as one would get from combining the populations into a single group: the groups are essentially identical and fully interact so are effectively a single group.

Additionally, the m-matrix is symmetric: if people in group a have a proportional rate of interaction with group b, then the group b has the same proportional rate of interaction with group a.

β matrix: how much people go into public

So I then make some bold assumptions:

A given β value is equally dependent on the behavior of the susceptible population and the relevant infectious population: the net β is proportional to the product of how much the S population interacts and how much the infectious population interacts.
The susceptible population S behaves similarly to the pre-symptomatic population P, which doesn’t generally realize it’s been infected.

So with these assumptions, each group’s S population interacts proportional to √β₀, and each group’s P population interacts also proportional to √β₀, so the rate at which P infects S within the group is the product of these terms, which is β₀.

There may also be differences in how infectious members of P are, but for purposes here, I fold that effect into the proximity component, which is essentially indistinguishable. So if people start wearing masks, for example, they are less likely to infect others, which is the equivalent of less proximity. But I assume that masks reduce infections comparably if worn by either party, which may be a poor assumption.

β matrix: likelihood encounter leads to infection

I still need to deal with other levels of infection within each group. I assume that these populations also interact with the outside world proportional to that group’s √β₀, and thus to get the correct β within the group, they have an infectiousness proportional to β₁/β₀ for I1. Since infection is the product of three parts: how much the susceptible population is in public (proportional to √β₀), how much the infectious population is in public (proportional to √β₀), and how infectious the population is (for I1, proportional to β₁/β₀), I get the correct β = β₁ multiplying these three terms. This is a big assumption, but I need to assume something to simplify the β matrix.

A similar approach is taken for groups I2 and I3: just replace in the above paragraph β₁ with β₂ or β₃.

I did it this way because by using a common “in public” weighting factor for each stage of the infection, I was able to simplify the calculation of how likely members of different groups would encounter each other, without doing separate calculations for each subset based on various β values (β₀, β₁, β₂, β₃). This was the first approach I tried and the math worked out relatively simply.

β matrix: probability of intergroup interaction

This is one part of β, but I cannot use these values directly, because I need to consider the rates of interaction between groups. There’s two aspects of interaction: more interaction between groups a and b result in a greater likelihood that members of b infect members of a, and vice versa, but the intermixing of the populations also implies a reduces probability that members of a infect members of a, or members of b infect other members of b. This is required by my second principle, that if groups are identical, I should be able to merge them and get the same result.

So consider groups a and b. If a susceptible member of group a encounters someone, where will be a probability pab that the encountered person will be from the opposite group b, and a probability (1 — pab) that the person will be from the same group a. The probability pab is a combination of how many members are in each group, how likely each member of the groups are to be in public, and how much the two groups intermix.

Multiplying these three terms, I get the following:

pab =( mab Nb √β₀b) / ( Na √β₀a + mab Nb √β₀b )

This is the probability that someone encountered by group a will be in group b. There is no maa term here, because maa = 1 implicitly: people interact with members of their own group.

β matrix: combining components

Combing these components yields the following formula, for two groups a and b. It is readily extensible to more than two groups. My code implementation allows for an arbitrary number of groups, with different mixing scales mab for the different group pairs.

d Ea / d t =
Sa [ Na √β₀a ( β₀a Pa + β₁a I1a + β₂a I2a + β₃a I3a) +
mab Nb √β₀b ( β₀b Pb + β₁b I1b + β₂b I2b + β₃b I3b) ] /
[ Na √β₀a + mab Nb √β₀b ]

d Eb / d t =
Sa [ mab Na √β₀a ( β₀a Pa + β₁a I1a + β₂a I2a + β₃a I3a) +
Nb √β₀b ( β₀b Pb + β₁b I1b + β₂b I2b + β₃b I3b) ] /
[ mab Na √β₀a + Nb √β₀b ]

I define the Sa, Sb, Ea, Eb, Pa, Pb, I1a, I1b, I2a, I2b, I3a, I3b to be fractions of their respective groups, each able to vary from the interval [0, 1]. In contrast, Na and Nb are total populations in groups a and b.

Note I just use a single interaction strength mab, since I assume interaction strength is symmetric. I don’t have a separate mba.

It’s evident that these equations meet my criterion that if all of the β values in the two groups are the same, and if intermixing is 100% (mab = 1), these groups can be combined into a single group with population Na + Nb, where the population fractions Sa, Sb, Ea, Eb ,Pa, Pb, I1a, I1b, I2a, I2b, I3a, I3b are combined into S, E, P, I1, I2, I3 fractions representing the combined populations.

running the model

I assume the β values are the same for each group except for scale factors for each group which may change over time due to interventions, or perhaps seasonal variation, etc.\
I then assign m-matrix values to define the interactions between groups. These could also be variable over time, for example due to travel restrictions.
I assign populations to each group.
I set up an initial infection. This can be in one group or in multiple groups. It’s popular to start a simulation with a small fraction of E in one or more groups, representing an initial exposure, or I like to use the eigenmode of steady-state exponential growth, which results in an instant transition to exponential increase.
Each iteration, I calculate the rates of change for each parameter for each group, including interactions, using parameter values at the start of the interval. I then use this rate of change to calculate the average value of each parameter during the interval, assuming either linear or exponential growth depending on the relative rate of change. I iterate this four times, so the rate of change and the average value are self-consistent. Finally I calculate the values at the end of the interval. I then move to the next time interval. A time interval of 0.1 days seems to work well. Exponential extrapolation, where appropriate, reduces the need for a very fine time step.

next step: results

I show some modeling results in a follow-up article.