A toy model for charting Malaysian COVID-19 trajectory

The SIR (susceptible/infected/removed) model in epidemiology is first considered nearly 100 years ago by two Scottish doctors, William Kermark (1898-1970) and Anderson McKendrick (1876-1943) (see W. O. Kermack and A. G. McKendrick (1927) Contribution to the mathematical theory of epidemics, Proc. Roy. Soc. Lond A 115, 700-721). In 1927, they opened their paper with the following line:

One of the most striking features in the study of epidemics is the difficulty of finding a causal factor which appears to be adequate to account for the magnitude of the frequent epidemics of disease which visit almost every population...

When a group of infectious species, \(I\), is introduced into a community of susceptible species, \(S\) (which is not previously exposed to the RNA materials), the virus spreads from the infectious species to the susceptible species, thus diminishing the susceptible species and growing the size of the infectious species.

$$I + S \longrightarrow I' + S' + R' \to R$$

Fortunately, a fraction of the infectious species will eventually recover from the disease, when they develop natural immunity (or receive forced-immunity via vaccination, the vaccine for COVID-19 is probably 12-18 months away) for fend off the virus attack and become noninfectious.

ICU and death summary (last updated on April 2)

Unfortunately, death is inevitable for a small fraction (probably 200 ppm or less, depending on the overall preparedness and load handling ability of the healthcare infrastructure in the community) of the infectious group, and collectively, the recovered/dead species are denoted by \(R\), since they are “removed" from the infectious group.

The SIR model, following the notation of Harko et al. (2004), is a set of ordinary differential equations:

\begin{align} \frac{dx}{dt} &= -\beta x y, & x(0) = N_1\\ \frac{dy}{dt} &= \beta x y - \gamma y, & y(0) = N_2\\ \frac{dz}{dt} &= \gamma y, & z(0) = N_3 \end{align}

where the symbols \(x, y, z\) denote the concentration or dimensionless size of the various interest groups (namely, \(S, I, R\), respectively). Any chemical engineer will immediately recognise that the problem considered is very similar to chemical kinetics, in which the transmission kinetics is first order in \(I\) and \(S\), while the recovery/removal kinetics is first order in \(I\).

I am a medical statistician. I've studied this stuff at university, done data analysis for decades, written several NHS guidelines (including one for an infectious disease), and taught it to health professionals. That's why you don't see me making any coronavirus forecasts.

— Robert Grant (@robertstats) April 2, 2020

It is quite easy to solve the SIR model with a suitable numerical tool, for example, with \(N_0 = 3.237\times 10^7, \beta = \frac{1}{3}, R_0 = 3, \gamma = \beta/R_0\) and with \(N_2 = k_2 N_2'/N_0, N_3 = k_3 N_3'/N_0\), where \(k_2 \approx 100, k_3 \approx 50\) are the multiplying factors introduced to “adjust" the reported statistics (since the reported numbers are only a subset of the actual size of the two groups of interest) by the Ministry of Health, e.g. as of March 29 noon, the reported numbers are:

$$N_3' = 388 + 34 = 422, N_2' = 2470 - N_3' = 2048$$

and with \(t = 0\) set to coincide with March 18, 2020 (the day of commencement of Movement Control Order), we predict that the size of the infectious group will reach its peak on Day 31 (that is, around April 18).

However, since the real values of \(R_0, \beta, \gamma, k_2, k_3\) are not exactly known (remember the tweet by Robert Grant), this peak value is just an educated guess based on the values initiated in the numerical simulation. The good news is, however, the gradient of the infectious trajectory is beginning to show some signs of slowing down.

For more toy models, see the excellent Python simulation by 3Brown1Blue, and the animated log-log plot by Aatish Bhatia and Henry Reich.