Infectious Diseases | Modelling Complex Systems | COVID-19 | Epidemiology |SIR Model | Virus Spread

#### Disclaimer

This series of articles was written in early 2020 and therefore most of the views and assumptions are outdated given our current state of knowledge.

## 1. Introduction

The SIR model has been historically used in epidemiology literature as a theoretical basis for understanding the establishment and spread of infectious diseases. It will, however, require a major overhaul if we were to apply it to the COVID-19 pandemic we are experiencing at the time of writing of this article.

In Part 4 of this series, we will attempt to leverage our previous work in Part 2 and Part 3 and the knowledge we have at our disposal from online data on the CDC website, Wikipedia, and medical journals to create a rigorous mathematical model.

Once the model is well defined, we will run simulations of different settings to try and understand the dynamics of the spread in various scenarios and discuss any noteworthy observations.

## 2. SEIR-COVID Model

### 2.1 Model Description

The SEIR-COVID model presented in this article is founded on the following ideas:

- The model is an extended version of the SEIR compartmental model
- The new model will accommodate additional groups as well as time delays

The seven groups described hereafter represent the different states in which an individual in the population can be found in at any point in time.

#### 2.1.1 Susceptible (S)

An individual starts in the Susceptible state, where he/she has not yet been exposed to the virus and has not developed immunity against it.

#### 2.1.2 Immune/Recovered (R)

People who have contracted the disease and who have not died will be moved to the Immune/Recovered state after a certain “recovery” period.

Given how important the role which **age**, and possibly other factors, play in the percentage of immune individuals in a given population, individuals can either start in the Susceptible state or in the Immune/Recovered.

#### 2.1.3 Exposed (E)

Exposed individuals are transferred from the Susceptible group to the Exposed group. The time spent in this state is known as the **incubation period**. During this stage of the disease, the person is not contagious and shows no symptoms.

#### 2.1.4 **Infected **(I)

When the incubation period is over, individuals will move to the Infected group during which the virus becomes transmissible by close contact with susceptible individuals.

#### 2.1.5 **Self-Isolation **(L)

The SEIR-COVID model supports a Self-Isolation state by allowing a percentage “s” of the infected population to move into that state.

Individuals in Self-Isolation are, by definition, not in close contact with other individuals and would not play any role in the generation of new cases.

For a closer representation of real life scenarios, we assume that there is a time delay associated with this stage, i.e. between the end of the incubation period and the start of self-isolation. This can be explained by:

- The failure of an individual to get tested as soon as symptoms are observed
- The inherent delay in the processing of the PCR test and communication back to the infected individual

#### 2.1.6 **Hospitalization** (H)

A percentage of infected people will eventually require hospitalization and would be moved to the Hospitalized group.

Given the criticality of available **ICU beds** and **ventilators** in the context of COVID-19, it is assumed, by definition, that only people who require intensive care will belong to this group.

#### 2.1.7 **Dead **(D)

The outcome of an infection is either death or recovery. It is assumed that a percentage of the population will die after contracting the disease.

Patients requiring intensive care but unable to have access to it will be assumed to die the next day and are moved to the Dead state.

In this model, we have assumed that people would not die in the Infected or Self-Isolated state if the Hospitalization group is activated (by having the model parameter “h” strictly greater than zero).

### 2.2 Transition between Groups

The below diagram shows the transition paths between the different groups. The parameters of the model: r, s, h, d, dh, as well as the delays will be examined closely in the coming sections. For the purpose of this diagram, suffice to know that these parameters represent the percentage of people moving from one group to the next.

At the moment of writing, it is not yet confirmed whether immunity can be lost over time. We will assume at this stage that recovered individuals will retain their immunity for the duration of the pandemic, which we will limit to one season or one life cycle.

### 2.2 Definitions and Notations

The following variables will be used in the subsequent sections:

- N = total number of individuals in the population
- S = number of individuals in the susceptible group
- E = number of exposed individuals
- I = number of infected individuals
- L = number of individuals in self-isolation
- R = number of recovered individuals
- D = number of individuals who died

Each individual will belong to one and only one of the seven groups at any point in time.

We also define the following model parameters:

- = contact rate between individuals of the same cluster
- = fraction of confirmed infected individuals going into self-isolation
- = fraction of infected individuals requiring intensive care
- = fraction of people who have died from the disease
- = fraction of people who have died while being hospitalized

We define the time delays as follows:

- = incubation period
- = period between end of incubation and onset of symptoms. This delay will determine the start of self-isolation
- = recovery period from end of incubation
- = period between end of incubation and death
- = period from end of incubation till hospital admission
- = period from start of hospitalization till recovery
- = period from start of hospitalization till death
- = the maximum number of available ICU beds for a population of 100,000.

### 2.3 Model Dynamics

The following diagram shows the transition of an exposed individual between the different states.

The set of differential equations governing the dynamics of this model are presented in this section.

First, we start with the rate of change of the **Susceptible** group which can be defined as:

Where:

This expression requires some explanation:

- The number of new infections is given to be directly proportional to the ratio of infected individuals in the population, . The proportionality ratio is also known as the
**Contact Rate**. - The contact rate is a time-dependant variable. This is important for the model to remain stable given the arbitrary nature of the constant . Also, as the number of immune individuals increases, the probability of an infected individual to pass the virus on to a susceptible individual would necessarily become smaller.

To simplify the remaining equations, we define as:

The number of people currently **Exposed **can be defined as:

- The number of people who just got infected minus the number of people who just completed their incubation period:

The rate of change in the **Infected** group can be defined as:

- The average number of exposed people who have completed their incubation period,
- Minus the number of people who have either started self-isolation, were admitted to hospitals, or have either died or recovered.

Applying the same concepts to the **self-isolation** group:

Following the same path, the rate of change in the number of **Hospitalized** people can be described by the following equation:

The rate of change in the **Recovered** group can be given by:

Similarly, for the rate of change in the **Death **group:

### 2.4 Population Changes

Given the relatively short duration of the epidemic, vitality and mortality rates will be ignored. Therefore, during this period, it will be assumed that the population size will not change and the total number of people in all groups will always be equal to the initial population size. This is represented by the following equations:

And:

### 2.5 Clusters

It might prove beneficial to look at the population not as a single large set of closely interacting individuals but more or less as groups of smaller sizes that we will call clusters. Basically, clusters can be defined as groups of individuals in closer contact with each other as compared to two individuals belonging to two different clusters. This structuring has been fully detailed in Part 2.

A cluster can be defined as a group of individuals categorized according to the seven groups we have defined for COVID-SEIR:

where and = the number of clusters in the system

A **lockdown **in this context would mean a severed interaction between two clusters. **Social Distancing** can be defined as a smaller value of the contact rate between individuals belonging to the same cluster.

The effect of clusters can be shown through the addition of a cross-contamination component which we will call .

The rate of change of S can thus be viewed as a result of inter- and intra-cluster contamination.

The term defines the cross-contamination from neighboring/interacting clusters on cluster . It can be written as:

Where:

- represents the
**Inter-cluster Contact Rate**between clusters and and is a random number between 0 and 1. The higher this number, the stronger the chance of cross-contamination between the two clusters. - is defined as the Lockdown Function. If a lockdown is enforced, there is a high chance that this interaction is severed and then , otherwise .

### 2.5 Cap on Ventilators and ICU Beds

The availability of ICU beds in conjunction with ventilators plays a major role in the fatality ratio of COVID-19. This will be factored in the model as follows:

- It is assumed that within a population of 100,000 individuals, a certain number of ICU beds with ventilators will be made available for COVID-19 patients. This is represented in our model by the constant . In some countries, this number can be as high as 38 for every 100,000 individuals.
- At time , the number of patients in ICU is and the number of expected patients in ICU can be given by .
- If , we assume that individuals will die.

Let be the expected number of new patients requiring hospitalization at . It is given by:

We define the fatality ratio as:

Where is the Heaviside step function.

Although the expected number of new hospital admissions would be , in practice, this number will be smaller due to the cap imposed on H(t):

The differential equations governing D(t), R(t), and H(t) would then look as follows:

### 2.6 Shorter Recovery Period

The recovery period (or duration of viral carriage) might prove to be important in the context of the spread of COVID-19. A shorter recovery has the potential of reducing hospital admission and thus preventing the NHS from being overwhelmed. It also means that less people would be likely to die as they would not move to the serious or critical stages of the disease which often happens on days 10-12 and above (see CDC Website on clinical progression).

In the context of this model, a shorter recovery period has an implication on the value of “h”, the percentage of infected individuals requiring hospitalization. We redefine the value of “h” as follows:

Where U(x) is the step function of Heaviside.

Although no treatment as yet has been widely adopted and approved by the World Health Organization, two protocols are being fiercely debated by the scientific community. One is Hydroxychloroquine+Azythromycin and another is Ivermectin+Doxycycline. Both of these drugs have shown a reduced duration of viral carriage.

### 2.7 Adding a Stochastic Element

The equations detailed in the previous section assume that the delays (admission to hospital, onset of symptoms, recovery, etc.) are fixed with zero variation. This is not entirely in agreement with the physical world.

To remedy this situation, we will assume that the delays experienced in each group follow an exponential distribution of averages where .

The modified equations take the following form:

Where the decay rate in the exposed group is:

The removal rate in the infected group is:

With:

Also, being the fraction of new patients requiring hospitalization and for which ICU care is not available.

While the removal rate in the self-isolation group is:

The removal rate in the hospitalized group is:

## 3 Model Simulation

As in Part 3 of this series, Python code was developed to generate numerical simulations of a physical system governed by the differential equations presented above. The Python code used to produce the simulations is available on this link.

As in the previous part, we define a baseline scenario and then vary the different parameters to produce variations for comparison and assessment.

### 3.1 Key Parameters

#### 3.1.1 Asymptomatic Cases

In this research paper, it was found that around 51% of people who tested positive in a specific setting (a cohort of sheltered homeless individuals in public facilities, asylum seekers, and employees in these facilities) did not have fever or respiratory difficulties. The section on Asymptomatic and Pre-Symptomatic Infection in the CDC website confirm this number.

#### 3.1.2 Immune Cases

In the same research paper, it was concluded that around 7% of residents and employees in housing facilities for homeless people tested positive for SARS-COV-2.

The cruise ship setting provides similar conclusions. The below table contains summary statistics from four cruise ships which had infected individuals onboard. The records included in this table were selected due to the high level of testing performed on the crew and passengers.

Ship | People | Tested | Cases | Deaths |
---|---|---|---|---|

Diamond Princess | 3711 | 3618 | 712 | 14 |

Celebrity Apex | 1407 | 1444 | 224 | 0 |

Horizon | 250 | 250 | 125 | 0 |

Greg Mortimer | 217 | 217 | 128 | 1 |

Totals | 5585 | NA | 1189 (21%) | 15 (1.26%) |

The significance of this information lies in the fact that residents had access to shared kitchens, toilets, bathrooms, and open spaces. This data could put a cap on the ratio of possible infections in an unmitigated scenario. There does not seem to be any data (aside from age) to explain the lower ratio of infections.

In reality, the median age of the crew and passengers on cruise ships could be lower than what would be observed in a city or suburb. This number was 35 for the Ruby Princess.

In our simulations, we have therefore used a conservative estimate of 40% for the population initially immune rather than 79%.

#### 3.1.3 Ratio of Self-Isolation

For a better representation of the real world, we assume that not all infected individuals move straight into self-isolation as soon as the incubation period is over. Following are the arguments for this hypothesis:

- We have assumed that the onset of symptoms occurs after a certain number of days beyond the incubation period. During this stage of the infection, the infected individual shows no symptoms but is otherwise contagious.
- Evidence shows that a number of people who have contracted the disease are asymptomatic and therefore are contagious but not aware of it.
- For different reasons, infected people have not strictly adhered to self-isolation rules and therefore continued to contribute to the spread of the disease.

It is difficult to quantify the value of “s”. From 3.1.1, we know it has to be less than 50%. Assuming a quarter of infected individuals fail to self-isolate, we will set “s” to a value of 40%.

#### 3.1.4 Cap on ICU Beds

In any real life scenario, the number of ICU beds per capita is limited. As this is critical to our model, we have again used a conservative estimate which puts this number at 35 per 100,000 inhabitants. This was based on data from Wikipedia.

A critical nuance is also the number and distribution of ventilators. For example, Russia has only 8 ICU beds vs 27 ventilators per 100,000 inhabitants. In this case, 35 looks to be a valid number as well.

Given that ICU beds would need to be made available for patients with other diseases, we will assume in the following simulations that only 24/100,000 beds are to be used for COVID-19 patients.

#### 3.1.5 Fatality Ratio

The Fatality Ratio parameters ‘d’ and ‘d_{h}‘ represent the percentage of infected individuals who die at home and in hospital respectively. Looking at the Worldometers website and in particular at those countries which

- Had the highest number of tests performed per 1M inhabitants, and
- Which are now late in pandemic, we observe the following:

Country | Total Cases | Deaths | Tests /1M | Population |
---|---|---|---|---|

Italy | 247,832 | 35,146 | 113,698 | 60,453,777 |

Spain | 335,602 | 28,445 | 142,834 | 46,756,402 |

UK | 303,952 | 46,193 | 239,707 | 67,917,167 |

Germany | 211,077 | 9,226 | 95,530 | 83,807,610 |

Belgium | 69,402 | 9,845 | 144,107 | 11,594,051 |

Totals | 1,167,865 | 128,855 | 147,145 | 270,529,007 |

The average fatality ratio for this group is:

- Dead / Total cases = 11%
- Dead / 1M = 476 (or 0.047%)

We will make the following assumptions in our simulations:

- Infected individuals who advance to critical stages in their disease are immediately hospitalized and therefore d is set to 0
- ICU patients form 6% (see CDC website) of the infected population. Therefore,
- Of all patients in ICU, around 50% are assumed to die. This is a worst case scenario as ICU care has improved dramatically in the last months. For our simulations, we will set .

### 3.2 Baseline

The model parameters used as our baseline are outlined below.

```
params = {
'r' : 2.3, # Contact rate
'd' : 0.00, # Death rate
's' : 0.4, # Self-isolation rate
'h' : 0.06, # Hospitalization rate
'dh' : 0.5, # Death rate in hospital
'tauE' : 5, # Incubation period
'tauR' : 20, # Recovery period
'tauD' : 15, # Time to die
'tauS' : 2, # Onset of symptoms
'tauH' : 12, # Hospital admission
'tauDH' : 8, # Hospital discharge on death
'tauRH' : 8, # Hospital discharge on recovery
'maxEpochs' : 2000,
'E0' : 1, # Initially Exposed
'S0' : 6000, # Initially Susceptible
'M0' : 4000, # Initially Immune
'H_max' : 24, # Max number of ICU beds per 100,000
'growthFunction' : 'Exp', # Linear, Gompertz, Constant, Exp
'lockdownActiveCases': -1, # Active cases before social distancing applies
'liftLockdownActiveCases': -1 # Active cases when social distancing stops applying
}
```

Metric | Value | Percentage |
---|---|---|

Peak (days) | 45 | NA |

Duration (days) | 240 | NA |

Susceptible (K) | 0.52 | 5.2% |

Immune (K) | 8.94 | 89.4% |

Death Rate | 543 | 5.4% |

### 3.3 Impact of Self-Isolation

This simulation shows the evolution of the spread in a scenario with no self-isolation.

```
params_adj['s'] = 0
```

Metric | Value | Baseline |
---|---|---|

Peak (days) | 35 | 45 |

Duration (days) | 250 | 240 |

Susceptible (K) | 1% | 5.2% |

Immune (K) | 93.5% | 89.4% |

Death Rate | 5.5% | 5.4% |

Observations:

- There does not seem to be any notable difference between the baseline and the scenario with no self-isolation. This could be explained by the fact that self-isolation is only applied to a small percentage (40%) of the entire infected population

### 3.4 Social Distancing

Social distancing means a smaller contact rate (of the order of 0.2 rather than 2.3) applied once the active cases reach a certain (arbitrary) limit, in this case 100.

The following parameters apply:

```
params_adj['lockdownActiveCases'] = 100
```

Once the number of infections (individuals in ICU plus individuals in self-isolation) reach the threshold, the lockdown will kick in and the Contact rate would change from to 0.2. In this simulation, the lockdown is not lifted until the virus is eradicated.

Metric | Value | Baseline |
---|---|---|

Peak (days) | 20 | 45 |

Duration (days) | 175 | 240 |

Susceptible (K) | 55% | 5.2% |

Immune (K) | 44.6% | 89.4% |

Death Rate | 0.04% | 5.4% |

Observations:

- The death rate is tiny compared to the baseline
- The higher ratio of susceptible individuals (55%) persisting after the disease has been eradicated implies the possibility of a second wave as we shall see in the next section.

### 3.5 Cyclical Lockdowns

In this simulation, lockdowns are imposed and lifted according to two parameters:

```
params_adj['lockdownActiveCases'] = 150
params_adj['liftLockdownActiveCases'] = 150
```

Once the number of infections (individuals in ICU plus individuals in self-isolation) reaches the given value for ‘lockdownActiveCases’, the lockdown will kick in and the Contact rate would be changed from to 0.2. If the number of identifiable infections decreases until it reaches ‘liftLockdownActiveCases’, the lockdown will be lifted.

Metric | Value | Baseline |
---|---|---|

Peak (days) | Multiple | 45 |

Duration (days) | 700 | 240 |

Susceptible (K) | 15.5% | 5.2% |

Immune (K) | 81% | 89.4% |

Death Rate | 3.5% | 5.4% |

Observations:

- The pandemic duration more than doubled, implying a heavier toll on the economy and well-being of the populace
- The death rate decreases by 35% compared to the baseline

### 3.6 Unlimited ICU Beds

In this simulation, we explore the evolution of the system when the cap on the number of available ICU beds is lifted.

```
params_adj['H_max'] = 9999
```

Metric | Value | Baseline |
---|---|---|

Peak (days) | 40 | 45 |

Duration (days) | 240 | 240 |

Susceptible (K) | 5.2% | 5.2% |

Immune (K) | 92% | 89.4% |

Death Rate | 2.8% | 5.4% |

Observations:

- Compared to the baseline scenario, the death rate almost halved due to the availability of more ICU beds.

### 3.7 Shorter Recovery Period

In this experiment, the impact of a shorter recovery period (decreased to 10 days) is explored. As mentioned earlier, we assume that the availability of a treatment would remove the need for any future hospitalization of infected individuals.

```
params_adj['tauR'] = 10
```

Metric | Value | Baseline |
---|---|---|

Peak (days) | 40 | 45 |

Duration (days) | 140 | 240 |

Susceptible (K) | 6% | 5.2% |

Immune (K) | 94% | 89.4% |

Death Rate | Nil | 5.4% |

Observations:

- The duration of the pandemic is shorter by 42%
- There is no possibility of a second wave as virtually all individuals have attained immunity
- The death rate is nil but thats more a result of the design than an outcome of the current settings

### 3.8 Higher Contact

The disease has been described as highly infectious but what does that mean for its evolution and spread? It is interesting to observe the changes in a setting where the contact rate has been elevated from 2.3 to 10, which is the objective of this simulation.

```
params_adj['r'] = 10
```

Metric | Value | Baseline |
---|---|---|

Peak (days) | 20 | 45 |

Duration (days) | 200 | 240 |

Susceptible (K) | 1% | 5.2% |

Immune (K) | 93.1% | 89.4% |

Death Rate | 5.9% | 5.4% |

Observations:

- The infection would run through the population in a slightly shorter time frame
- The death rate is a shade higher (9%) as compared to the baseline

### 3.9 Impact of Clusters

This numerical experiment examines the impact of clustering individuals in smaller groups. The figure on the left shows which clusters are in contact while the heatmap on the right shows the strength of the interaction between two clusters.

The model parameters for this scenario are:

```
nbrOfClusters = 10
clusterInteraction = 0.2
clusterInteractionVar = 0.1
```

We will consider three cases: a baseline, single lockdown, and cyclical lockdown.

#### 3.9.1 Baseline with Clustering Enabled

Observations:

- The dynamics of the spread of the disease and how it ends are not very different from what we have seen for the single group cases
- All 10 clusters have been infected at the early stages of the spread
- The spikes in the Day-on-Day graph indicate the date at which a new outbreak has been recorded in a new cluster

#### 3.9.2 Lockdown with Clustering

In this scenario, we observe the system while a lockdown is applied after a certain threshold in identifiable infections has been reached.

```
params_adj['lockdownActiveCases'] = 50
```

Obervations:

- No significant deviations observed between this setting and the previous one in section 3.4
- The lockdown, under these settings (contact rate between clusters follows a normal distribution with average 0.2 and standard deviation of 0.1) did not produce significantly different results. All the clusters were subsequently infected.
- The death rate is significantly lower than the baseline

#### 3.9.3 Cyclical Lockdown with Clustering

This numerical experiment allows for the lifting of the lockdown when the number of identifiable infections reaches a certain lower limit, in this case set to 50.

```
params_adj['lockdownActiveCases'] = 50
params_adj['liftLockdownActiveCases'] = 50
```

Observations:

- The duration of the pandemic is significantly higher (around 8 times)
- All the clusters were infected
- The fatality ratio is halfway between the baseline and the single lockdown
- The dynamics anticipate several subsequent waves until 79% of the population is immune.

## 4 Conclusion

Experience has shown time and time again that using models to predict or even speculate about the future is a very risky endeavor. We believe that this holds true for the COVID SIR model presented here as well.

The complexity of the ecosystem formed by the virus, its host, and the transmission channels is far too complicated to be accurately represented in a set of simple differential equations with a little more than a few parameters. The unpredictability and sometimes irrationality of human behavior takes the problem to a different level.

What modelling does provide however are powerful insights into the dynamics of an otherwise opaque and complex dynamic system. They would allow us to separate some wheat from the chaff as the latter would necessarily produce outcomes that do not agree with real world observations.

Finally, when the pandemic is over (or stabilizes into an equilibrium state), the model can be updated or fine-tuned and its parameters could be determined with a higher level of precision. This may then allow scientists a more transparent view on this subject.