Modelling the Spread of COVID-19: Part 3: Simulation Results

1. Simulation Exercise

In Part 2 of this series, we have spent some time defining the software model that we will be using to simulate the spread and impact of the coronavirus throughout a susceptible population.

In this article, we present the details of the model parameters. We will also combine these parameters to reproduce scenarios (Eradication, Suppression, Herd Immunity, Test-Diagnose-Treat, etc.) from the real world. We then observe the results of these simulations.

The following metrics can be quite important and need to be analyzed as functions of the model parameters:

• Peak number of infections
• Duration of the epidemic
• Fatality rate as a function of age, number of ICU beds, and recovery period
• Herd immunity
• Truth around the second wave
• The case for the lockdowns
• Impact of Social Distancing

The Python code used to produce the simulations is available on this link.

2. Model Parameters

The parameters used as input to the model are described below. The majority of the values used in the simulations have been collected from the CDC website.

• populationSize : 50. This is the average size of a single cluster. For the purpose of this simulation, we consider clusters to represent groups of people who are in close contact long enough for virus transmission to occur. These can be households or small size working places such as restaurants or pharmacies.
• epochs : 1000. This is the maximum number of epochs (or days) for which the simulation will run.
• numberOfClusters : 200. Total number of clusters in the system.
• clusterSocialDistanceType : ‘Low’. The starting point in the simulation for Social Distancing. It is considered that social distancing within closely connected small groups such as families cannot be observed for practical reasons.
• interClusterSocialDistanceType : ‘Low’. This parameter represents the social distancing factor as observed between clusters. This parameter changes to Medium or High depending on the scenario settings, such as when the Lockdown is activated.
• contactRate : 1.3. The contact rate represents an arbitrary number greater than one. This is not to be confused with the Basic Reproduction Rate or R₀.
• daysToRecover : 20. This number represents the average number of days till recovery. It follows a lognormal distribution with mean 0 and var 0.4.
• daysToRecoverWithTreatment : 0. This number represents the average number of days till recovery if a treatment is being administered to the patient. It follows a lognormal distribution with mean 0 and var 0.4. If this parameter is not 0, it will used for non-asymptomatic individuals.
• selfIsolation : True. This parameter enables or disables the self-isolation feature in the software model. Self-isolation happens on the onset of symptoms. If the value is set to True, on every cycle, the infected person has a 90% chance of going into self-isolation provided that he is not asymptomatic.
• incubationPeriod : 5. The incubation period is the delay between exposure to the virus and onset of symptoms. It follows a lognormal distribution centered around 5 with mean 0 and variance 0.4. During the incubation period, it is considered that the exposed individual is not contagious.
• percAsymptomatic : 15. The percentage of asymptomatic people is controlled via this parameter. It only applies to individuals below 40 years old.
• lockdownActive : False. This parameter simulates the effect of a lockdown. If its set to True, interactions between clusters have a 90% chance of being severed allowing the prevention of cross-contamination between two clusters.
• liftLockdownCases : 50. Number of active cases before a lockdown is lifted. This is more of an arbitrary number as there is no science to back up any such decisions.
• activateLockdownCases : 50. Number of active cases before a lockdown is initiated. This is more of an arbitrary number as there is no science to back up any such decisions.
• interactions. The matrix of interactions indicates which clusters participate in cross-contamination. It is generated using a uniform distribution where Pr(0)=4.8/5 and Pr(1)=0.2/5.
• interactionsWeights. A matrix defining the interaction weights between two clusters. This represents the percentage of individuals from Cluster A coming into close contact with individuals from cluster B.
• childrenImmune : True. At the time of writing, immunity among children and young adults appears to be well established. If this parameter is set to True, individuals younger than 20 years old have been considered immune, individuals between 20 and 30 years old are 50% of the time immune, and individuals above 30 years of age are considered susceptible.
• maxHospitalCapacity : 35. This the maximum number of ICU beds per 100,000 inhabitants. The value of 35 chosen for the simulations represents one of the highest in the world.
• daysToICUAdmission : 8. Average number of days between onset of symptoms and admission to ICU.

The death of an infected individual could occur under the following circumstances:

• If the individual’s age is above 80, chance of dying from the disease is 15%.
• If the individual’s age is between 70-79, the fatality rate is 8%
• If the individual’s age is between 60-69, the fatality rate is 4%.
• If the individual’s age is between 50-59, the fatality rate is 2%.
• If the individual has been infected for more than daysToICUAdmission and is above 50 years of age, the probability of requiring intensive care in a hospital is 4% for age bracket 50-59, 8% for 60-69, and 16% if above 80 years.
• If the maximum ICU capacity indicated by maxHospitalCapacity has been reached and the individual is not able to be admitted into a hospital, its assumed that the individual will die immediately.

3. Simulation Scenarios

3.1 Generic Setup

The following parameters have been setup to form the baseline for comparison:

nbrOfClusters = 100
clusterInteraction = 0.2
clusterInteractionVar = 0.1
averageAge = 35

params = {
    'populationSize' : 100,
    'epochs' : 1000,
    'numberOfClusters' : nbrOfClusters,
    'clusterSocialDistanceType' : 'Low',
    'interClusterSocialDistanceType' : 'Low',
    'contactRate' : 2.3,
    'daysToRecover' : 20,
    'daysToRecoverWithTreatment' : -1,
    'selfIsolation' : True,
    'incubationPeriod' : 5,
    'percAsymptomatic' : 50,
    'lockdownActive' : False,
    'liftLockdownCases' : 20,
    'activateLockdownCases' : 50,
    'interactions' : None,
    'interactionsWeights' : None,
    'childrenImmune' : True,
    'MaxHospitalCapacity' : 24,
    'daysToICUAdmission' : 12
}

The heatmaps depicted below show the interactions between the different clusters. Clusters A and B are connected if the value of Interactions(A, B) = 1. The number of connections can be controlled via a binomial distribution where Pr(0)=4.8/5 and Pr(1)=0.2/5.

The age distribution can be seen in the following histogram. The latter is skewed to the left to simulate a younger population which is more common in developing countries or larger cities.

3.2 Baseline

The baseline simulation can be characterized as follows:

• The population has some level of immunity (37%) in children and young adults but is otherwise susceptible
• No policies (lockdowns, social distancing, treatment, etc.) except for self-isolation are imposed

The purpose is to:

• Understand the dynamics of the spread in the simplest setup
• Form a baseline against which variations of this scenario can be compared

After running the simulation with the above parameters, we observe the following results:

All four graphs show one or more properties of the system as measured against time, the latter being measured in days:

• The graph on the top left corner shows the evolution of the infectious population with time.
• The one on the top right corner shows the immune and susceptible populations
• The diagram on the bottom left corner shows the number of clusters infected with at least one individual.
• Finally, the bottom right graph shows the number of deaths as well as the total years lost. This is defined as 80 minus age of the individual at the time of death. The number 80 was chosen to represent the lifetime expectancy in most countries.

Observations:

• The pandemic ended only when 100% of the population was immune

3.3 Impact of Self-isolation

The purpose of this simulation is to:

• Understand the impact of a policy with no self-isolation.

The parameters for this model are the same as the baseline except for self-isolation which in this case is set to False:

params_adj['selfIsolation'] = False

The following diagram shows the result of the simulation:

The following observations can be made:

Conclusion:

• There does not seem to be a difference in the overall results with or without self-isolation.
• This can be explained by the fact that virus transmission continues to occur in asymptomatic cases (50% of the overall infections), infected individuals who do not self-isolate (25%), and infected individuals between incubation and onset of symptoms

3.4 Higher Contact Rate

In this scenario, we try to:

• Understand the impact of a higher contact rate, which can be interpreted in this context as a higher contagiousness of the disease or a denser population.

This can be controlled by increasing the value of the contact rate from 1.3 to 4.

params_adj['contactRate'] = 4

The following diagram shows the result of the simulation:

We have the following observations:

Conclusion:

• The overall situation does not appear to change when the contact rate is higher (or much higher) than 2.3.
• The peak of the pandemic occurs earlier as compared to the baseline and the overall duration is slightly shorter
• The death rate decreases by 11.1%

3.5 Shorter Recovery Period

One of the proposed treatments of COVID-19 is the drug combination of Hydroxychloroquine + Azythromycine. It has the twin advantage of lowering the viral carriage duration and avoiding the need for hospitalization. The purpose of this simulation is to:

• Test the impact of a shorter recovery period

This scenario can be tested by decreasing the recovery period:

params_adj['daysToRecoverWithTreatment'] = 7

The following diagram shows the result of the simulation:

The following observations can be made:

Conclusion:

• The number of deaths significantly lowers by 56%

3.6 Eradication Approach

The Eradication approach consists of imposing a hard lockdown until the disease is eradicated.

The purpose of this scenario is to:

• Understand the results of this approach in the framework that has thus far been used (peak infections, death rates, duration of the epidemic)

This is achieved by having active cases drop to zero before lifting the lockdown. During the lockdown phase, the social distancing factor observed within a cluster is Medium (High is not practical as it means almost no contact between individuals within a cluster).

params_adj['lockdownActive'] = True
params_adj['clusterSocialDistanceType'] = 'Medium'
params_adj['activateLockdownCases'] = 150
params_adj['liftLockdownCases'] = -1

The below diagram shows the result of the simulation:

We have the following observations:

Conclusions:

• The overall duration of the pandemic is double that of the baseline
• The death rate decreases by 39.9%

3.7 Suppression Approach

The Suppression approach consists of keeping the disease under control by imposing and lifting the lockdown when necessary. The numbers set to control the lockdowns (lower bound for activation = 50 cases, higher bounds for lifting the lockdown = 50) are arbitrary numbers as there is no scientific theory or study to backup the case.

The purpose of this scenario is to:

• Understand the impact of a cyclical lockdown
• The nature of the second wave

params_adj['lockdownActive'] = True
params_adj['activateLockdownCases'] = 150
params_adj['clusterSocialDistanceType'] = 'Medium'
params_adj['liftLockdownCases'] = 150

The below diagram shows the result of the simulation:

We have the following observations:

Conclusions:

• The duration of the epidemic is significantly increased
• The simulation outcome shows a sequence of waves with a downward trend in its peaks
• The death rate is smaller by 17.6%

3.8 Medium Social Distancing

The final scenario presented here is one where only medium social distancing is imposed, i.e. without any lockdowns, much like the Swedish response.

The purpose of this simulation is to:

• Allows us to understand whether a lockdown can be avoided through social distancing alone.

The parameters used are:

params_adj['clusterSocialDistanceType'] = 'Medium'

The below diagram shows the result of the simulation:

We have the following observations:

Conclusions:

• The end result is not very different from the baseline.

3.9 Infinite ICU Capacity

This setting is designed to observe the impact of ICU beds on the overall evolution of the system.

The parameters used are:

params_adj['MaxHospitalCapacity'] = 9999

The below diagram shows the result of the simulation:

The following observations can be made:

Conclusions:

• The end result is not very different from the baseline.
• The death rate is 15.6% lower

4. Can Models be Used for Making Predictions?

Although this exercise cannot be considered a “model” in the traditional sense, the outcome of the numerical experiments and those generated from a model should be the same given that they both operate on the same underlying principles, the only difference being in how those concepts are embodied (software code vs set of differential equations).

Going back to the question, the short answer is No, as has been brutally demonstrated time and over again. The long answer is elaborated in the following thoughts.

4.1 Models as Embodiment of Ideas

Models are the embodiment of ideas that we believe to be true about the physical system we are trying to analyze.

For example, in the context of COVID-19, if we believe that children do not spread the disease and are immune to it, then we build the model to produce exactly that effect. The results will only mirror the rules and constraints that had produced them. Models are never expected to uncover the laws governing the dynamics of the system or propose alternatives to existing ones.

4.2 Over-simplification of the Reality

The model that we saw uses 15 adjustable parameters and around 10 hard-coded ones. This is a ridiculously small number as compared to the degrees of freedom that are presented in a real world physical system.

To realize a practical model that can be run on a computer in a reasonable amount of time, a trade-off in complexity (and therefore precision) must be made in the design.

In addition to higher computation times, another potential issue can arise with a higher number of “degrees of freedom” in a model. The model can become extra sensitive to small changes in these parameters due to the non-linearity of the internal forces governing it. It would also be very challenging, given the number of options one might end up with, to pick combination that would mimic the real world.

4.3 Modeling Human Behavior

Human behavior is only rational and predictable in hindsight but it can always be irrational and unpredictable in unforeseen situations. This adds an extra level of complexity when modeling interactions either between two human beings, or between humans and their environment.

People adapt their behavior to changing circumstances in different ways, ways that are sometimes difficult to comprehend. This can, for example, be seen in the people’s attitude towards self-isolation; although marketed heavily as a life-saver, it has been confirmed that it has not been followed properly.

4.4 Reality as a Constantly Moving Target

The world we live in today has never been so interconnected, dynamic, and fragile. Behavior and thought are heavily modulated by daily news broadcasts and social media. It is more the realm of chaos theory rather than deterministic, even stochastic, models.

Capturing the nature of reality with constant parameters is practical but not ideal. The internal laws governing the system change as time passes, or as a function of other forces, internal or external to the system.

When a model is conceived, it is implicitly assumed that the rules governing its dynamics will remain constant for some time. The reality however could be different and by the time the model is ready, the rules may have changed drastically.

5. What do Mathematical Models Really Offer?

Given all the inherent “deficiencies” in model-making and the apparent inability to derive precise predictions from models, one might ask the question of why do scientists spend such an enormous time and effort trying to build them?

Here are some thoughts that might help answer that question.

5.1 Testing Our Knowledge

Once a model is built, it could be tested with data and statistics collected from the real world. If the outcome agrees with observations, it would then be safe to assume that the underlying theories are reasonably valid.

In a way, modeling can be used to eliminate invalid theories or assumptions.

5.2 Understanding the Dynamics of the System

This is vital for any real understanding of the behavior of the virus and its impact on the population as we have seen when running the different scenarios. More crucial still, are the insights derived for policy making.

A good understanding of any hidden interactions would help formulate a better response to any future pandemic (its probably too late for this one).

Whats Next?

In the following articles, we will attempt to produce the differential equations governing the model we have simulated and compare the results with these simulations.