Second wave simulations

Introduction

This page accompanies my article on COVID-19 second wave published in The Conversation, by providing simulation details.

The aim of this model is not to predict the type, size and timing of a potential “second wave” of the COVID-19 epidemic, but rather to provide an illustration of potential scenarios. As such, it is kept very simple.

Model

The model used to produce the graphs is a basic SIR model with a time-dependent rate of spread. The population of $N$ individuals is split into four subpopulations, susceptible $S$ , infected $I$ , recovered $R$ and dead $D$ . Initially, all are susceptible except a small proportion of individuals who are infected. There is no attempt to model latency or different stages in the disease.

The model is described in terms of four differential equations:

$\begin{array}{l*{10}{c}c} \frac{dS}{dt}=-\beta(t) S \frac{I}{N} \\[10pt] \frac{dI}{dt}=\beta(t) S \frac{I}{N} - g I \\[10pt] \frac{dR}{dt}=g (1-p) I \\[10pt] \frac{dD}{dt}=g p I\end{array}$

where $\beta(t)$ is the time-dependent rate of spread, $g$ is the rate of removal, and $p$ is the proportion of deaths. Thus, $\tau=1/g$ is the length of the infectious period.

The time-dependent rate of spread is given by a combination of double logistic functions forming a pulse,

$f(t)=\frac{1}{1+\exp\left(-r (t-t_0+\Delta/2)\right)}\times \frac{1}{1+ \exp\left(r (t-t_0-\Delta/2)\right)}.$

where $r$ determines how fast the rate changed in response to lockdown or release of lockdown, $t_0$ is the centre of the pulse, and $\Delta$ is the width of the pulse.

The corresponding value of the current reproduction number is given by

$R_e=\frac{\beta(t)}{\tau}.$

The asymptomatic spread was modelled by assuming that only one in eight cases is symptomatic, i.e. the reporting rate is 12.5%.

For comparison with data, the number of infected individuals on any given day, $I(t)$ was converted into reports by firstly dividing by the asymptomatic rate (8) and subsequently dividing by the infectious period.

Parameter values

The infectious period $\tau=8$ is a value widely published for COVID-19 data. The death probability $p=0.019$ was estimated by comparing the output with the UK death data.

The values for the time-dependent rate are as follows:

	$r$ 1/day	$t_0$ days	$\Delta$ days	Max $R_e$	Min $R_e$	Comments
1	0.25	-21	100	2.7	0.8	First wave until lockdown
2	0.25	500	800	0.8	1	Rebound until $R_e=1$
3	0.25	500	800	0.8	1.2	Rebound until $R_e=1.2$
4	0.25	(120, 220, 320, 420)	50	0.8	1.2	Series of pulses with interventions
5	0.25	250	50	0.8	4	Autumn wave with a mutated strain

The initial value of $R_e$ is quite well established in the literature; it also agreed well with the initial growth in the UK , as shown in the graph below. The value of $R_e=0.8$ during the lockdown has also been attested by various methods for the UK data in late April and May. The increase to 1.2 was chosen to illustrate the danger of getting wrong the post-lockdown spread of the virus. The uncertainty of $\pm 20$ % is routinely reported for the values of $R_e$ .

The length of the pulses (50 days) reflect a combination of 7-14 days from the start of the increase in the rate to the detection of the increase by measuring $R_e$ , an implementation period for the future lockdown, and the period it takes for the lockdown measures to damp down the value of $R_e$ . 50 days just simply “looked good” by comparing the period with the “original” March-May 2020 lockdown.

The increase of $R_e$ for the autumn wave simulation was chosen to represent a high bound for the current estimate for SARS-CoV-2 ( $R_e=4$ ).

The probability that an infected individual dies, $p$ is equal to 0.019, i.e. 2%. This is broadly in agreement with data.

The initial number of all infected cases is $I(0)=1200$ and the simulations start on 5th March 2020. Other initial values were $S(0)=70,000,000-I(0)$ , $R(0)=0$ , $D(0)=0$ . Assuming only one-eighth of cases are symptomatic, this corresponds to 150 initially infected symptomatic individuals and to 19 reported cases (dividing further by the infectious period to convert infections to reports). The UK reported 48 cases on 5th March 2020.

Results

Equations were solved in R using the odeSolve package. The parameters were chosen to represent UK data until 26th May 2020, as shown below, although no formal fitting procedure was carried out. As I intended to keep the article in The Conversation as generic as possible, I did not compare the UK data with the model output, but the comparison is shown below. Lines are model input (rate) and output (reported cases and deaths), points represent data (Johns Hopkins University, collected on 28th May 2020).

The effective reproduction number, the number of reported cases and the total number of deaths as predicted by the model (solid line) compared to data for the UK (blue dots). The number of cases is estimated using the reporting rate of 12.5% and by dividing the number of infected individuals by the infectious period. No formal fitting was carried out. The period of the main lockdown is marked in grey; the start is 23rd March and the end is 23rd May as for England (Scotland, Northern Ireland and Wales implemented the lockdown for longer).

The number of reported cases and the daily number of deaths comparison. The number of cases is estimated using the reporting rate of 12.5% and by dividing the number of infected individuals by the infectious period. No formal fitting was carried out. Red line corresponds to model and data in complete agreement.

Caveats

The model is too simple to be able to predict the value of $\beta(t)$ at any given period as well as to fully account for the impact of the post-lockdown increase in $R_e$ . It cannot also be used to predict the number of cases or deaths in the future.

In the article, I did not want to go into details about simulations of deaths and did not show the comparison.

Thus, all the results in the article in The Conversation are for illustration purpose only.

Added 18th June

I have created a simulation tool that allows the user to explore the simulation model described above. It uses a slightly different parameterization for the start and end of the waves but should be self-explanatory. It will download the current UK data from the Johns Hopkins web site, so at some point, the parameters will stop working. Please use the menu to select the options and press Update button to run. The full-page version is available here.