SIR: Factors Influencing Spread¶

REPLICABILITY

Randomness is in an important factor in RKnot’s approach to simulation (and frankly in real-word viral transmission), so the results of the sims below will not be repeatable with each iteration. The below examples are meant to show general differences based on state; further analysis should run the same simulation multiple times to see the mean impact.

Base US Simulation¶

To explore the various concepts of RKnot and viral spread, we’ll use a simulation design based on CDC Best Planning Scenario guidelines for COVID-19 characteristics including:

\(R_0\) 2.5
IFR for each of 4 age groups
- 0-19: 0.003%
- 20-49: 0.02%
- 50-69: 0.5%
- 70+: 5.4%

Other assumptions:

Population of \(10,000^1\)
Initial Infected of 2
Duration of Infectiousness 14 days\(^2\)
Duration of Immunity 365 days
Density of ~1 subject per location (dlevel='med')

\(^1\)proportionately split among the 4 age groups to match US census data.

\(^2\)equal likelihood of transmission on any day (i.e. no viral load curve)

US Census Data

Natural¶

1. Equal¶

The first simulation makes the most homogeneous assumptions.

No group is restricted in terms of movement.
All dots are able to interact with one another.
All dots are susceptible at initiation.
All dots are equally likely to move to any dot on the grid (mover='equal')

The basic layout is below. These parameters can also be imported from rknot.sims.baseus for convenienve.

group1 = dict(
    name='0-19',
    n=2700,
    n_inf=0,
    ifr=0.00003,
    mover='equal',
)
group2 = dict(
    name='20-49',
    n=4100,
    n_inf=1,
    ifr=0.0002,
    mover='equal',
)
group3 = dict(
    name='50-69',
    n=2300,
    n_inf=1,
    ifr=0.005,
    mover='equal',
)
group4 = dict(
    name='70+',
    n=900,
    n_inf=0,
    ifr=0.054,
    mover='equal',
)
groups = [group1, group2, group3, group4]
params = {'dlevel': 'med', 'Ro':2.5, 'days': 365, 'imndur': 365, 'infdur': 14}

We instantiate a new sim by passing groups and params. We can also flag details to get some information about the Sim structure.

from rknot import Sim, Chart
sim = Sim(groups=groups, details=True, **params)
sim.run()

---------------------------------------------------------------------------------
|                                  SIM DETAILS                                  |
|-------------------------------------------------------------------------------|
|           Boundary|      [ 1 45  1 45]|          Locations|              2,025|
|-------------------|-------------------|-------------------|-------------------|
|         Population|             10,000|            Density|               4.94|
|-------------------|-------------------|-------------------|-------------------|
|       Contact Rate|               4.94|                   |                   |
|-------------------|-------------------|-------------------|-------------------|

Running the animation will result in a video as per below:

chart = Chart(sim, use_init_func=True)
chart.animate.to_html5_video()

Note embedded videos are used for convenience purposes. Given the random processes involved, running the same code will produce slightly different results each time.

Results:

Peak	37%
HIT	67%
Total	87%
Fatalities	0.55%
% > 70	42%
IFR	0.63%
Days to Peak	72

In this scenario, RKnot fairly closely replicates the curve of a standard SIR model, which expects HIT of 60% for \(R_0\) of 2.5 (HIT = 1 - 1 / \(R_0\)).

Variation from the standard SIR model will always result given:

In the simulation, movement and transmission are stochastic processes.
This sim does not have an entirely homogeneous population, with different IFRs and varying numbers of contacts between subjects.

2. Local¶

In remaining simulations, we begin to introduce ever increasing homogeneity.

Our first change is to adjust the subjects mover functions to local. The local mover has a strong bias towards locations only in its immediate vicinity, which is a better approximation of real world processes (though certainly not a perfect one).

We will also extend the simulation an extra year.

group1['mover'] = 'local'
group2['mover'] = 'local'
group3['mover'] = 'local'
group4['mover'] = 'local'
groups = [group1, group2, group3, group4]
params['days'] = int(365*2)

sim = Sim(groups=groups, **params)
sim.run()

chart = Chart(sim, use_init_func=True)
chart.animate.to_html5_video()

Results:

Peak	4%
HIT	56%
Total	69%
Fatalities	0.42%
% > 70	36%
IFR	0.60%
Days to Peak	355

Compared to Example 1, we can see that restricting movement has a major impact on spread. The curve is flattened and extended with total infected, peak, HIT, and fatalities all reduce.

But the virus is never completely eradicated and instead moves in a progressive wave across the grid space. Note that the infection and fatalitiy levels are somewhat artificially reduced as the virus still has a large susceptible population in the upper portion of the grid that it has not yet reached.

Local movement does not satisfy the “well-mixed” condition of the SIR model. Particularly interesting is that the initial conditions of this scenario *DO* very closely resemble a SIR model. So, an observer measuring \(R_0\) in a local environment might see spread consistent with a well-mixed population but the population may not be well-mixed at a larger scale.

4. Pre-Existing Immunity¶

In this scenario, we adjust the susceptibility factor for just two groups by relatively small amounts as follows:

20-49: 80%
50-69: 65%

This means, in the inverse, that 10% and 25% of the subjects in the respective groups are already immune to the virus (whether through pre-existing T-cell immunity, anti-bodies, or otherwise).

The older group is assumed to have a lower susceptibility factor as it is more likely that older people will have had more exposure to similar viruses over their lifetime.

T-cell immunity to sars-cov-2 remains a controversial subject, but many studies have found prevalance of T-cells between 20% - 50% in people unexposed to sars-cov-2. It is suggested that exposure to “common cold” coronaviruses (or more dangerous ones) may convey this immunity.

group2['susf'] = .8
group3['susf'] = .65
groups = [group1, group2, group3, group4]

sim = Sim(groups=groups, **params)
sim.run()

chart = Chart(sim, use_init_func=True)
chart.animate.to_html5_video()

Results:

Peak	19%
HIT	44%
Total	63%
Fatalities	0.35%
% > 70	29%
IFR	0.56%
Days to Peak	93

Relative to Example 3 above, pre-existing immunity would reduce all aspects of the spread curve signifcantly.

Note the reduction in fatalities results even though the most susceptible group is assumed to NOT have pre-existing immunity.

5. Events¶

Certainly, the vast majority of people in the US do not move in such pre-defined ways as set out by the mover function. In reality, people tend to move with a local bias with a small number of interactions, supplemented by larger movements to locations with a large number of interactions in a small amount of time.

In RKnot, we can simulate this with Event objects. And we will incorporate a number of them in this simulation.

First, we will reset our group parameters by importing from sims.baseus.

Then we will instantiate a handful of events recurring periodically over the duration of the sim.

from rknot.sims.baseus import params, groups

from rknot.events import Event

school1 = Event(
    name='school1', xy=[25,42], start_tick=2,
    groups=[0], capacity=10, recurring=2
)
school2 = Event(
    name='school2', xy=[78,82], start_tick=3,
    groups=[0], capacity=10, recurring=2
)
school3 = Event(
    name='school3', xy=[92,32], start_tick=4,
    groups=[0], capacity=7, recurring=2
)
game1 = Event(
    name='game1', xy=[50,84], start_tick=6,
    groups=[0,1,2,3], capacity=100, recurring=14
)
game2 = Event(
    name='game2', xy=[45,52], start_tick=5,
    groups=[0,1,2], capacity=76, recurring=14
)
game3 = Event(
    name='game3', xy=[12,87], start_tick=1,
    groups=[0,1,2], capacity=56, recurring=14
)
game4 = Event(
    name='game4', xy=[52,98], start_tick=3,
    groups=[1,2], capacity=113, recurring=28
)
concert1 = Event(
    name='concert1', xy=[20,20], start_tick=7,
    groups=[0,1], capacity=50, recurring=14
)
concert2 = Event(
    name='concert2', xy=[91,92], start_tick=28,
    groups=[1], capacity=50, recurring=14
)
concert3 = Event(
    name='concert3', xy=[62,38], start_tick=21,
    groups=[2,3], capacity=25, recurring=14
)
concert4 = Event(
    name='concert4', xy=[38,42], start_tick=14,
    groups=[1,2], capacity=50, recurring=14
)
bar1 = Event(
    name='bar1', xy=[17,24], start_tick=4,
    groups=[1], capacity=5, recurring=3
)
bar2 = Event(
    name='bar2', xy=[87,13], start_tick=5,
    groups=[1], capacity=5, recurring=4
)
bar3 = Event(
    name='bar3', xy=[52,89], start_tick=6,
    groups=[1,2], capacity=5, recurring=3
)
bar4 = Event(
    name='bar4', xy=[16,27], start_tick=7,
    groups=[1,2,3], capacity=4, recurring=7
)
bar5 = Event(
    name='bar5', xy=[89,46], start_tick=6,
    groups=[1,2], capacity=7, recurring=7
)
church = Event(
    name='church', xy=[2,91], start_tick=7,
    groups=[2,3], capacity=20, recurring=7
)

events = [
    school1, school2, school3, game1, game2, game3, game4,
    concert1, concert2, concert3, concert4,
    bar1, bar2, bar3, bar4, bar5, church
]

sim = Sim(groups=groups, events=events, **params)
sim.run()

chart = Chart(sim, use_init_func=True)
chart.animate.to_html5_video()

Results:

Peak	29%
HIT	56%
Total	82%
Fatalities	0.50%
% > 70	39%
IFR	0.61%
Days to Peak	77

We’ve attempted to tailor the event setup so that the infection curve resembles that of Example 1. The theory is that the \(R_0\) measured in the early stages of a pandemic should be reflective of actual contacts, not the theoretical contacts of the model.

This scenario is not perfectly substitutable with Example 1, however.

One measure we can use to compare scenarios is the total number of interactions. We find that under the Event and Gated approaches it takes more contacts to result in the same level of spread.

Avg Number of Contacts per Subject:

Example 1: 50.3

Example 2: 52.3

Example 3: 53.4

Example 5: 65.5

Example 6: 75.5

Average of first 100 days across 5 sims for each scenario.

Other issues with substitutability may exist and are being explored. SIR models determine R0 in the idealized environment of Example 1 and so may not be suitable for customized environments such as this one.

6. Gates¶

We can further improve the real world relevance of interactions by introducing gates. Subjects are not always freely able to interact with all other people in a population. Often their movement is restricted to within certain areas. Furthermore, other people’s access into those areas is restricted.

Elderly people living in retirement homes or assisted living centers is an example. To simulate this, we will split group4 into two separate groups.

group4a
- population of 600 (2/3s of group4)
- IFR of 4.2%
- move freely throughout the entire grid as previously
group4b
- population of 300 (1/3rd of group4)
- IFR of 7.8%
- movement restricted to 6x6 box

We have also adjusted the IFR on the basis that group4b is likely older and also probably more frail than group4a. IFRs approximate those found here.

In addition, we will add an event specifically for the new group inside the gate.

group4a = dict(
    name='70+',
    n=600,
    n_inf=0,
    ifr=0.042,
    mover='local',
)
group4b = dict(
    name='70+G',
    n=300,
    n_inf=0,
    ifr=0.0683,
    mover='local',
    box=[1,6,1,6],
    box_is_gated=True,
)
groups = [group1, group2, group3, group4a, group4b]

church2 = Event(
    name='church2', xy=[2,3], start_tick=7,
    groups=[4], capacity=5, recurring=7
)

events_gated = [
    school1, school2, school3, game1, game2, game3, game4,
    concert1, concert2, concert3, concert4,
    bar1, bar2, bar3, bar4, bar5,
    church, church2,
]

Now, such elderly populations are entirely sealed of from the rest of the world. In fact, they are often visited by family members or friends. We can mimick this with the use of a Travel object.

In this sim, at least one person will enter into the group4b gate for a day only. And this will repeat every day of the sim.

from rknot.events import Travel

visit = Travel(
    name='visit', xy=[1,1], start_tick=3,
    groups=[1,2], capacity=1, duration=1, recurring=1
)
events.append(visit)

sim = Sim(groups=groups, events=events, **params)
sim.run()

chart = Chart(sim, use_init_func=True)
chart.animate.to_html5_video()

Results:

Peak	30%
HIT	62%
Total	83%
Fatalities	0.52%
% > 70	40%
IFR	0.62%
Days to Peak	66

Again, the gated structure is intended to mimic Example 1, Example 3, and Example 5.

The use for Events and Gates will become clear when we explore the impact of Policy Decisions.

7. Pre-Immunity with Events and Gates¶

Now we can see how pre-immunity might impact viral spread in a population with more heterogeneous interactions.

group2['susf'] = .8
group3['susf'] = .65
groups = [group1, group2, group3, group4a, group4b]

sim = Sim(groups=groups, events=events, **params)
sim.run()

chart = Chart(sim, use_init_func=True)
chart.animate.to_html5_video()

Results:

Peak	20%
HIT	42%
Total	62%
Fatalities	0.47%
% > 70	37%
IFR	0.76%
Days to Peak	86

Similar to Example 4, pre-existing immunity flattens the curve somewhat and results in a modest decrease in fatalities (even with an abnormally high IFR in this case).

10. Isolation¶

TBD

Policy¶

With a more realistic model of subject interaction, we can begin to experiment with the impact of different policy measures.

RKnot can simulate policy measures via Restriction, SocialDistancing, and Quarantine objects. Further details here.

The simulations are based on this scenario and so the impact of the policy measures contemplated should be considered relative to that scenario.

The structure can be imported as follows:

from rknot.sims.baseus import params, events_gated, groups_gated

1. Restrict Large Gatherings¶

We’ll start by simply restricting large gathers, which for this sim is any event with 10+ capacity (0.1% of the entire population). The restrictions will last for 120 days.

When considering capacity, remember that a 100,000-seat stadium in a 10 million person catchment represents 1% of the population.

We assume this policy are implemented on day 30, when the population finally realizes there is a pandemic and government has had time to implement prevention measures.

The restriction will last for 120 days.

from rknot.events import Restriction

large_gatherings = Restriction(
    name='large', start_tick=30, duration=120, criteria={'capacity': 10}
)
events_w_res = events_gated + [large_gatherings]

sim = Sim(groups=groups_gated, events=events_w_res, details=True, **params)
sim.run()

chart = Chart(sim, use_init_func=True)
chart.animate.to_html5_video()

Results:

Peak	15%
HIT	28%
Total	73%
Fatalities	0.45%
% > 70	38%
IFR	0.62%
Days to Peak	69

We see that the curve is flattened significantly during the restriction period. Infections have a much lower peak, BUT also a much longer tail (evident by the still high relatively high level of total infections).

This results as the event restriction is lifted, which leads to a slower rate of decline of the virus.

3. Restrict Elderly Visits¶

We can restrict events by name by passing name key to criteria. We can do this to restrict the travel event to the 70+G area.

The policy measure is implemented on day 30 and maintained for another 120 days. There will be no other restrictions.

no_visits = Restriction(
    name='no_visits', start_tick=30,
    duration=120, criteria={'name': 'visit'}
)
events_w_res = events_gated + [no_visits]

sim = Sim(groups=groups_gated, events=events_w_res, **params)
sim.run()

chart = Chart(sim, use_init_func=True)
chart.animate.to_html5_video()

Results:

Peak	27%
HIT	54%
Total	80%
Fatalities	0.29%
% > 70	17%
IFR	0.36%
Days to Peak	77

By simply quarantining the elderly, we see the most dramatic reduction in fatalities of any scenario thus far. An outbreak NEVER occurs in the 70+G gated area, even AFTER the policy restrictions are lifted, because the virus has already been eradicated by that time.

We can see here that allowing a high level of infection among the least susceptible has resulted in a low level of fatalities among the most susceptible.

This is a controversial approach and does have difficult ethical implications, but its power to reduce death cannot be ignored.

5. Quarantine¶

We can even mimick the impact of quarantines via the Quarantine object.

Here we show the impact of a 30-day quarantine for all groups in the Sim. We include a restriction on visits to the elderly during the quarantine.

from rknot.events import Quarantine

quarantine = Quarantine(
    name='all', start_tick=30,
    groups=[0,1,2,3,4], duration=30
)
no_visits = Restriction(
    name='no_visits', start_tick=30,
    duration=30, criteria={'name': 'visit'}
)

events_w_res = events_gated + [quarantine, no_visits]

sim = Sim(groups=groups_gated, events=events_w_res, **params)
sim.run()

chart = Chart(sim, use_init_func=True)
chart.animate.to_html5_video()

Results:

Peak	1%
HIT	1%
Total	1%
Fatalities	0.00%
% > 70	0%
IFR	0.00%
Days to Peak	32

The pandemic never materializes in this simulation. The quarantine eradicates the virus swiftly upon implementation.

This is an interesting result. This obvioulsy did not occur in many places in the real world that took this approach. There are several possible reasons worth considering:

Random Variation: if you run this simulation multiple times, you will note some sims where spread does occur.
Scale: a larger simulation would increase the possibility that a very small number of subjects can continue to pass around the virus while it is muted in the broader population
Adherence: real-world adherence to the implemented policies was much lower than this simulation suggests.
Inconsistency: some areas implemented strict quarantines while others did not, and there was mixing among those populations.

6. Quarantine with Adherence¶

We can investigate the impact of low adherence to quarantine requirements by setting the adherence property of the Quarantine object.

from rknot.events import Quarantine

quarantine = Quarantine(
    name='all', start_tick=30,
    groups=[0,1,2,3,4], duration=30
)
no_visits = Restriction(
    name='no_visits', start_tick=30,
    duration=30, criteria={'name': 'visit'}
)

events_w_res = events_gated + [quarantine, no_visits]

sim = Sim(groups=groups_gated, events=events_w_res, **params)
sim.run()

chart = Chart(sim, use_init_func=True)
chart.animate.to_html5_video()

[16]:

chart_params

[16]:

{'use_init_func': True,
 'show_intro': True,
 'dotsize': 0.1,
 'h_base': 10,
 'interval': 100.0}

With just 10% of the population ignoring the quarantine requirements, the virus spreads almost as though there was no quarantine at all.

The virus barely survives for the first 6-months of the outbreak, then as per other scenarios, once restrictions are lifted an outbreak occurs. Still, when the outbreak does occur, it is characterized by one of the lowest peaks and HITs in our analysis (due in part to the pre-immmunity of some of the groups).

And it achieves the lowest fatality rate of the group, mainly by restricting access to the elderly population for the duration of the pandemic and ensuring an outbreak never occurs in that region.

SIR: Factors Influencing Spread¶

Base US Simulation¶

Natural¶

1. Equal¶

2. Local¶

3. Social¶

4. Pre-Existing Immunity¶

5. Events¶

6. Gates¶

7. Pre-Immunity with Events and Gates¶

8. Self Aware Social Distancing¶

9. Self Aware Social Distancing + Pre-Immunity¶

10. Isolation¶

Policy¶

1. Restrict Large Gatherings¶

2. Social Distancing¶

3. Restrict Elderly Visits¶

4. Social Distancing and Restrict Elderly Visits and Large Gatherings¶

5. Quarantine¶

6. Quarantine with Adherence¶