Modeling the COVID-19 Outbreak with J | DK’S ABODE

🚨

the banner says: trans rights!

Foreword: I am not an epidemiologist and I don’t claim to be. I am, however, a math major and an academic and I feel like I’ve done my due dilligence in reporting this accurately and correctly. The full code used is at the end of the post, if you’d rather not sift through data.

If you’ve been watching the news, I assume you’ve heard about the COVID-19 pandemic. In light of the CDC releasing a model predicting 62% of Americans will become infected1, I became interested in how they came to those numbers. In fact, how do we come to any conclusions as to the spread of disease at all?

There are plenty of well established models for disease spread. The simplest model is the SIR model5, but it comes with a few caveats. The first caveat is that it assumes that once a disease runs its course in a single person, that person will now be permanently immune to the disease. The second caveat is that the model assumes that all infectious people immediately begin to show symptoms. Neither of these assumptions are accurate for COVID-19: the disease has a 14% reinfection rate6, and it takes 10 to 14 days for symptoms to appear.

The model that has been used for the current COVID-19 pandemic is the SEIRS model*,4,5. With the SEIR model, we have 4 states that any particular person can be in: “Susceptible”, “Exposed”, “Infectious”, and “Recovered”7.

A person is susceptible if they are at risk for catching a disease.

A person is exposed if they have caught the disease or are carrying the disease without any symptoms. This happens during the incubation period.

A person is infectious if they have now caught the disease and are experiencing symptoms. This is when people go to the hospital, require intensive care, etc.

A person is recovered if they have gotten over the disease and now have some semblance of immunity.

We have 4 variables we can tweak: \(\beta, \sigma, \gamma, \xi \). They are all rates over a generic unit of time.

\(\beta\) is the rate at which susceptible people become exposed.

\(\sigma\) is the rate at which exposed people become infectious.

\(\gamma\) is the rate at which infectious people recover.

\(\xi\) is the rate at which recovered people lose immunity and become susceptible.

Used with attribution from https://www.idmod.org/docs/hiv/_images/SEIR-SEIRS.png. Copyright 2019, Intellectual Ventures Management, LLC (IVM). All rights reserved.

These states and variables all have differential equations which correlate them. There are still few notes to be made, though. Firstly, this model doesn’t represent mortality in the population due to the disease. Because of that, the basic reproductive rate4,7 (that is, the amount of people one person will infect, on average) of the disease is quite simple to calculate: \(R_0 = \frac{\beta}{\gamma}\).

We’re going to complicate things a little bit, though. Just looking at disease numbers is boring and unimpactful! Instead of analyzing the differential equations, we’re going to take a look at an actual discrete time simulation which uses those 4 variables. But first, we have to figure out what those variables should be.

Let’s take our “unit of time” to be one day. Here’s some of the information that we have for the virus:

We know that \(R_0\) is about 2.288 for COVID-19.

We know that it has an incubation period of about 10-ish days (2-14, as per the CDC).

We know that once symptoms emerge, they persist about two weeks.9

We know that 14%6 of infected people become re-infected.

We can calculate all the variables we need from these, starting with \(\sigma\). Each day, there’s a certain chance \(\sigma\) for COVID-19 to begin showing symptoms. Each day, the chance that it cumulatively hasn’t shown symptoms is \((1 - \sigma)^N\), where \(N\) is the number of days. This is a simple binomial distribution, and for it to begin showing symptoms after 10 days on average means that we must simply solve for \((1 - \sigma)^{10} = 0.5\). So, \(\sigma = 0.066967\)

Once symptoms emerge, they persist about two weeks. Let’s do the same process we did above: \((1 - \gamma)^{14} = 0.5\) means we have \(\gamma = 0.0483048\).

We can calculate \(\beta\) using \(R_0 = \frac{\beta}{\gamma}\) with the \(R_0\) and \(\gamma\) values established above8. \(2.28 = \frac{\beta}{0.0483048}\) means \(\beta = 0.110134944\).

All’s left is to calculate \(\xi\). I have no clue how we could do this, but let’s just say (as a completely random, reasonable guess) that it takes about 30 days to lose immunity to COVID-19. \((1 - \xi)^{30} = 0.5\) gives us \(\xi = 0.02284\).

Now that we have all our variables, let’s boot up our J session and start writing this simulation!

First, numeric definitions.

ppl =: 10 NB. How many people are in our simulation?

Beta =: 0.110134944 NB. The rate at which susceptible people become exposed.<br>Sigma =: 0.066967 NB. The rate at which exposed people become infectious.<br>Gamma =:...