Stochastic compartmental models of the COVID-19 pandemic must have temporally correlated uncertainties

1. Introduction

Quantitative models for forecasting the spread of communicable disease are a valuable tool for policymaking during pandemics, such as the ongoing COVID-19, and endemic diseases, such as gonorrhoea. Compartmental models are simple, widely used, quantitative tools for studying the spread of such infectious disease [1–6]. For instance, an appropriate compartmental model for the COVID-19 pandemic is the susceptible-exposed-infected-removed (SEIR) model [7–15], depicted in figure 1a. As shown in figure 1c, if the model parameters are chosen properly, the deterministic SEIR model accurately forecasts the spread of the Omicron variant in the USA.

• Download figure
• Open in new tab
• Download PowerPoint

However, the model parameters, such as the infectivity rate $\lambda$, the incubation rate $\alpha$ and the curing rate $\gamma$, are a priori unknown. These parameters can be estimated by collecting and averaging data over the population. Such data are inevitably riddled with uncertainties that have motivated the development and use of stochastic compartmental models.

In stochastic compartmental models, one assumes that the model parameters are either random variables or stochastic processes [9,16]. The former assumes that parameter values are uncertain, but otherwise constant over time. In the present work, we model the parameters as stochastic processes, in order to account for random changes in parameters that are time-varying, and possibly correlated in time. Note that this time-varying uncertainty in model parameters is different from models with deterministic time-varying coefficients that account for the seasonal variability [17]. The most common approach in models with stochastic processes as parameters is to assume that the parameter comprises a constant mean perturbed with white noise [9,16,18,19]. For instance, the infectivity or contact rate is often modelled as $\lambda (t)=\bar{\lambda }+\sigma \xi (t)$, where $\bar{\lambda }$ is the mean, $\xi (t)$ is the standard white noise and $\sigma$ is a constant controlling the noise intensity [16,20–25]. Here, we first show that, although this assumption seems reasonable, it leads to systematic underestimation of the disease spread. For instance, figure 1d shows the SEIR prediction of the Omicron wave in the USA with the infectivity rate $\lambda$ modelled as white noise. The resulting stochastic SEIR model consistently underestimates the true cumulative number of COVID-19 cases. Making matters worse, as the noise intensity $\sigma$ increases, the white noise stochastic model further underestimates the number of COVID-19 cases. The mean trajectories of COVID-19 cases drop significantly compared with the predictions of the deterministic model. For $\sigma =0.6\bar{\lambda }$, the deterministic forecast is even outside the region where the 50% of trajectories around the mean trajectory lie (figure 1d).

As we discuss in §3, this behaviour is not specific to the SEIR model, but is also observed in SIR and SIS models. This indicates that the stochastic compartmental models with white noise have a fundamental flaw: they imply that greater uncertainty in the model parameters leads to a less severe spread of the disease. This model behaviour is certainly dubious since greater ignorance about the disease, e.g. its infectivity rate, does not in reality make a pandemic less severe. The counterintuitive and unrealistic implications of modelling parameter uncertainties with white noise had also been noted in oncology models, where increasing the noise intensity leads to tumour eradication, a behaviour that is not replicated in actual cases [26].

By contrast, as shown in figure 1e, modelling uncertainties by correlated noise alleviates such contradictory behaviour. With correlated noise, stochastic compartmental models exhibit the behaviour one would expect: their mean trajectories are in line with the deterministic model for all values of $\sigma$, and higher noise intensity leads to higher variance around the mean trajectory.

Here, for the first time, we present a principled modelling of uncertainties in the infectivity rate of epidemiological models. Starting from reasonable assumptions on the contact rate of each individual, and using the central limit theorem and Doob’s theorem on Gaussian Markov processes, we prove that the only admissible model for such uncertainties is the mean-reverting Ornstein–Uhlenbeck (OU) process. We demonstrate the implications of our results on two examples, an SIS model of sexually transmitted disease and SEIR model of COVID-19. For the SIS model, we derive a closed-form solution for the probability distribution of the infected population and show that the OU-based model always has a peak near the deterministic equilibrium. By contrast, modelling uncertainties with white noise erroneously predicts a less severe disease spread, or even its eradication, as the uncertainties grow. Similar observations are made for the SEIR model of COVID-19 using Monte Carlo simulations and the Johns Hopkins University database of COVID-19 cases.

We note that many epidemiological models have already incorporated stochasticity [27]. A straightforward approach is the stochastic perturbation of the model parameters [16,20–25]. A more principled approach is to model the number of susceptible and infected individuals as continuous-time Markov chain birth–death processes or as branching processes [28–34]. Under the usual assumption that their state variables are continuous, Markov chain models lead to a system of stochastic differential equations whose noise covariance matrix is different from the parametric perturbation models [28,35–37]. Nonetheless, these equations use white noise to account for uncertainties, overlooking the temporal correlations present in human social activities [38].

The OU process has also been proposed for perturbation of parameters in compartmental epidemiological models [28,39–41]. So far, the justification for the use of OU noise has not been based on first principles. Rather, OU has been chosen for its mean-reverting property; on average, OU noise stays closer to its mean value compared with the white noise. Thus, by choosing OU instead of white noise, unrealistic values of model parameters, though not eliminated, are less frequent.

The first contribution of the present work is the observation that compartmental models with white noise fluctuations in contact rate systematically underestimate the severity of disease spread in the population. Secondly, we justify the use of OU noise fluctuations from first principles. More specifically, we prove that the OU process is the only admissible noise that takes into account the temporal correlations present in the social behaviour of individuals. Finally, we show that the predictions of compartmental models with OU noise are drastically different from the predictions of white noise models; OU model predictions are always in agreement with the actual data, and do not underestimate the disease spread.

2. Stochastic modelling of the average contact rate

There are several parameters in compartmental models of infectious disease. For instance, in the SEIR model, these are the infectivity or average contact rate $\lambda$, the incubation rate $\alpha$ and the curing rate $\gamma$. The source of their uncertainty is the underlying dynamical interplay between biological, social and environmental factors. Incubation and curing rates depend mainly on the biology of the virus. However, the contact rate $\lambda$, defined as the average number of adequate contacts per individual per unit time ([42], Sec. 2.7.1), [43], depends on both biological factors, e.g. how easily the virus transmits, and social factors, e.g. how many individuals each person comes in adequate contact with. As a result, the contact rate is the main source of uncertainty in epidemiological models. Therefore, for the sake of simplicity, we assume that the average incubation and curing rates are deterministic constants and that the contact rate $\lambda$ is the only stochastic parameter. This is a common choice in the literature on stochastic compartmental models [10,16,44].

We denote the total number of adequate contacts of the $n$th individual, up to time $t$, by $C_n(t)$. A contact is adequate if it is sufficient for the transmission of the virus [45]. The contact rate of the $n$th individual is then defined as the rate of change of $C_n$, i.e. ${\lambda }_n=\mathrm{d}{C}_n/\mathrm{d}t$. Since $C_n(t)$ is not necessarily differentiable, it is more convenient to work with its increments, $\mathrm{\Delta }{C}_n(t)$, that is, the number of contacts of the $n$th individual over the time interval $[t,t+\mathrm{\Delta }t]$, where $\mathrm{\Delta }t$ is a small time increment. The contact rate is then approximated by ${\lambda }_n(t)\simeq \mathrm{\Delta }{C}_n(t)/\mathrm{\Delta }t$.

We treat the incremental contacts $\mathrm{\Delta }{C}_n(t)$ of each individual as a stochastic process satisfying the following assumptions:

(1)	There is no dependence between the incremental contacts of different individuals. In other words, we assume ${\{\mathrm{\Delta }{C}_n(t)\}}_{n=1}^N$ is a collection of $N$ independent random processes.
(2)	The number of contacts each individual makes in $[t,t+\mathrm{\Delta }t]$ does not depend on the history of their contacts in previous time instances; thus $\mathrm{\Delta }{C}_n(t)$ is a Markov process, where $N$ is the population size.
(3)	The mean and covariance of the incremental contacts satisfy $$\mathbb{E}[\mathrm{\Delta }{C}_n(t)]=\mu \mathrm{\Delta }t$$2.1a and $$\mathrm{Cov}[\mathrm{\Delta }{C}_n(t),\mathrm{\Delta }{C}_n(s)]=\frac{{\sigma }_0^2}{2\tau }\exp (-\frac{|t-s|}{\tau })\mathrm{\Delta }{t}^2,$$2.1b where $[t,t+\mathrm{\Delta }t]$ and $[s,s+\mathrm{\Delta }t]$ are disjoint time intervals, and $\tau$ is the correlation time.

Equation (2.1) models the assumption that the incremental contacts of every individual are proportional to the time interval $\mathrm{\Delta }t$. Equation (2.1a), where $\mu$ is a constant, models the assumption of no seasonal variability in contacts. The exponential decay in equation (2.1b) constitutes the simplest model for temporally correlated incremental contacts. The correlation time $\tau$ can be interpreted as a characteristic time scale for the social patterns of every individual (e.g. one day or one week). Note that, although the incremental constants $\mathrm{\Delta }{C}_n(t)$ have the same mean and covariance, we do not assume that they are identically distributed.

Assumptions (1)–(3) determine the properties of the contact rate ${\lambda }_n(t)$. In particular, ${\lambda }_n(t)$ inherits the Markovian property of $\mathrm{\Delta }{C}_n(t)$. Also, the mean value and covariance of ${\lambda }_n(t)$ are determined via equation (2.1),

and The exponentially decaying correlation of equation (2.2b) for the contact rate of each individual was also assumed in the agent-based model of Ariel & Louzun [39]. Having defined the contact rate ${\lambda }_n(t)$ for every individual, we now introduce the contact rate averaged over the whole population as $\lambda (t)=(1/N)\underset{n=1}{\overset{N}{\sum }}{\lambda }_n(t)$.

Now, we consider $K$ time instances $t_1<\cdots <t_k<\cdots <t_K$ and define the vectors ${\mathit{\lambda }}_n={[{\lambda }_n(t_1),\dots ,{\lambda }_n(t_k),\dots ,{\lambda }_n(t_K)]}^T$, which samples the stochastic process ${\lambda }_n(t)$ at the said time instances. The random vector ${\mathit{\lambda }}_n$ has mean value $\mathit{\mu }$ with ${\mu }_k=\mu$, and covariance matrix $\mathit{\Sigma }$ with ${\Sigma }_{k\ell }=\mathrm{Cov}[{\lambda }_n(t_k),{\lambda }_n(t_{\ell })]$, $k,\ell =1,\dots ,K$. For the average vector $\mathit{\lambda }=(1/N)\underset{n=1}{\overset{N}{\sum }}\lambda {\lambda }_n$, the multi-dimensional central limit theorem [46] implies

where $N$ is the total population, and ${\mathcal{N}}_K$ is the $K$-variate normal distribution. Note that $\mathit{\lambda }$ can be viewed as the discrete samples of the average contact rate $\lambda (t)$, so that $\mathit{\lambda }={[\lambda (t_1),\dots ,\lambda (t_k),\dots ,\lambda (t_K)]}^T$.

Thus, the stochastic process $\lambda (t)$ of the average contact rate has the following properties:

(i)	It is temporally homogeneous, since for every sequence $t_1<\cdots <t_k<\cdots <t_K$, all the $K$-variate joint distributions of $\lambda (t_1),\dots ,\lambda (t_K)$ are invariant to translations in time.
(ii)	It is a Markov process, since it is the average of independent Markov processes [47,48].
(iii)	For every pair of two time instances, $t$ and $s$, its values $\lambda (t)$ and $\lambda (s)$ follow a bivariate Gaussian distribution.

Doob’s theorem on Gaussian Markov processes [49] implies that the only two possible stochastic processes with properties (i)–(iii) are the uncorrelated white noise and the correlated OU noise. White noise is also ruled out since we have already proven by the central limit theorem that $\lambda (t)$ and $\lambda (s)$ are not independent. Therefore, the only choice for $\lambda (t)$ is to be an OU process, with mean value and autocovariance with $\sigma ={\sigma }_0/\sqrt{N}$.

It is well known that the OU process, with mean (2.4) and stationary autocovariance (2.5), is generated by the stochastic differential equation,

where $W(t)$ is the standard Wiener process [50]. We note that, as the correlation time $\tau$ tends to zero, the uncorrelated white noise is retrieved ([51], Sec. 6.6). In §3, we use this limiting relationship to compare results from stochastic compartmental models under white and OU noise.

For notational simplicity, we write $\lambda (t)=\bar{\lambda }+\sigma \xi (t)$, where $\bar{\lambda }=\mu$ is its mean value, $\sigma$ is its noise intensity and $\xi (t)$ is a standard OU process with zero mean and autocovariance given by equation (2.5) with $\sigma =1$.

3. Results

(a) Stochastic SIS models

Before discussing the COVID-19 data, we examine the simpler SIS model where we derive the probability distribution of the infected individuals in closed form. This allows us to highlight the differences between the white noise and OU noise with no ambiguity. The SIS model is governed by the equations

and with the susceptible (S) and infected (I) compartments where $N=S+I$ is the total population. We assume that the curing rate $\gamma$ is a deterministic constant and the contact rate $\lambda (t)=\bar{\lambda }+\sigma \xi (t)$ is a stochastic process. We consider two types of noise $\xi (t)$: white noise and the standard OU noise. It is customary to express the results in terms of the infected fraction $X=I/N$. Note that since the total population $N$ is conserved, we have $S=N(1-X)$. Throughout the paper, the noise intensity $\sigma$ is reported as a percentage of the mean value $\bar{\lambda }$, in order to measure the extent of variability in relation to the mean. The parameter $\sigma /\bar{\lambda }$ is known as the standard deviation to mean value ratio or coefficient of variance.

For stochastic SIS models, we are able to obtain the final distribution of the infected fraction $X$ of the population in closed form (see §S2 of the electronic supplementary material). As seen in figure 2, the theoretical final distributions are always in excellent agreement with Monte Carlo simulations.

• Download figure
• Open in new tab
• Download PowerPoint

The final distribution depends on two-dimensionless parameters (see §S2 of the electronic supplementary material), the deterministic basic reproduction number $R_0=\bar{\lambda }/\gamma$ and the relative variance ${\sigma }^2/\bar{\lambda }$ of the noise. Thus, we are able to produce the bifurcation diagram shown in figure 2, where the possible final forms of the distribution are depicted. For $R_0>1$, the deterministic SIS model ( $\sigma =0$), has one stable equilibrium at $X=(\bar{\lambda }-\gamma )/\bar{\lambda }$, and one unstable equilibrium at $X=0$. However, the stochastic SIS model with white noise exhibits a richer asymptotic behaviour that includes parameter regions of bistability and regions where $X=0$ is stable.

We first focus on $R_0=1.4$, which is the relevant basic reproduction number for gonorrhoea [45]. As shown in figure 2, by increasing the white noise intensity, the infected fraction distribution undergoes a transition. For low noise intensity, the distribution exhibits one mode (region I of figure 2), close to the stable deterministic equilibrium (figure 2a). As the noise intensity increases, the distribution mode shifts towards zero (figure 2b). For more noise intensity, the distribution exhibits its mode at zero (region III), see figure 2c. Increasing the white noise intensity further results in the eradication of the disease from the population (region IV), since the distribution of the infected fraction $X$ becomes a delta function at zero (figure 2d). We emphasize that this is an unrealistic model behaviour since it implies that greater ignorance about the contact rate $\lambda$ will lead to disease eradication.

For more contagious diseases with basic reproduction number higher that 2, the final distribution undergoes a different noise-induced transition. At first, as the noise intensity increases, the mode of the distribution (region I) shifts towards values that are larger than the deterministic equilibrium. For even higher noise intensity (figure 2e), an additional mode at zero appears (region II), whose magnitude increases by further increase in noise intensity. This increase eventually results in the eradication of the disease from the population, since the distribution becomes a delta function centred at zero (region IV).

In figure 2a–d, we also plot the distribution obtained under an OU perturbation of the contact rate in the SIS model. As correlation time of the OU noise, we choose $\tau =0.5{(\bar{\lambda }-\gamma )}^{-1}$, where ${(\bar{\lambda }-\gamma )}^{-1}$ is the characteristic time scale of the deterministic SIS model (see appendix A).

We observe that, compared with white noise perturbation, correlation in noise has the effect of suppressing transition away from the deterministic equilibrium. For noise levels as high as 80%, the mode of the distribution remains close to the stable deterministic equilibrium. For unrealistically large noise levels (such as $\sigma =1.2\bar{\lambda }$ of figure 2d), the distribution under OU perturbation also exhibits an additional mode at 0. However, most of distribution mass is still located around the deterministic equilibrium, whereas, for the same noise levels, the SIS model with white noise perturbation predicts the eradication of the disease from the population.

(b) Omicron variant in the USA

Next, we study the Omicron variant of COVID-19 in the USA, using the SEIR model,

with the susceptible (S), exposed (E), infected (I) and removed (R) compartments. It shares the same set of parameters with the SIS model as well as the additional parameter $\alpha$ denoting the average incubation rate. We note that more complex models which include a vaccinated compartment, e.g. the susceptible-vaccinated-exposed-infected-removed (SVEIR) model [52], can also be used. Following earlier studies [7–15], here we treat individuals that are effectively vaccinated as part of the removed compartment and therefore use the simpler SEIR model.

Our main COVID-19 results are shown in figure 1. First, we determine the parameters of the deterministic SEIR model by fitting its trajectory to the data from the cumulative COVID-19 cases in the USA [53] for the period of December 2021–April 2022, during the Omicron wave (figure 1b,c). The fitting procedure is described in appendix A. Then we consider the stochastic contact rate $\lambda (t)=\bar{\lambda }+\sigma \xi (t)$, which is perturbed around its deterministic value $\bar{\lambda }=1.54$ days ${\hspace{0.17em}}^{-1}$. For the perturbation $\xi (t)$, we consider both white and OU processes, and perform Monte Carlo simulations of the respective stochastic SEIR models (figure 1d,e). The noise intensity is varied between 10% and 60% of the average contact rate $\bar{\lambda }$. For the OU noise, we also have to choose the value of the correlation time $\tau$. Unfortunately, empirical data on the temporal correlation of contact rates are scarce [38]. For the results shown in figure 1e, we choose one week as the correlation time, motivated by the weekly pattern in human activity [54–57]. The results are similar for shorter and longer correlation times (see §S3 of the electronic supplementary material).

As shown in figure 1d, increasing the intensity of both white and OU noise results in an expected increase in variance of the forecasted Omicron cases. However, we observe that, in the white noise case, increasing the noise intensity results in the mean to underestimate the final number of cases, predicting significantly less spread of the disease in the population. For instance, when the noise intensity is at 60%, the forecasted cumulative number of COVID-19 cases is underestimated by approximately 20%. This result is unrealistic, since higher uncertainty in the contact rate does not in reality result in a less severe pandemic.

On the other hand, the results in figure 1e for the correlated OU noise are reliable. The mean trajectories stay close to the deterministic prediction and the data. As the noise intensity increases, the predictions attain slightly higher mean values of the final number of cases, compared with the value predicted by the deterministic model. For instance, at 60% noise intensity, the forecasted cumulative number of COVID-19 cases increases by approximately 4%.

The mechanism by which the SEIR model under white noise underestimates the size of the COVID-19 pandemic is better understood by inspecting the noise-induced qualitative changes in the distribution of the final number of cases; see figure 3. In absence of a closed-form solution, we use Monte Carlo simulations to determine the distribution corresponding to the stochastic SEIR model.

• Download figure
• Open in new tab
• Download PowerPoint

As shown in figure 3, for small noise intensities, the distributions of both white and OU stochastic models are unimodal and narrow, with their modes approximately coinciding with the deterministic equilibrium, marked by a dashed line. The distribution under OU noise is slightly less diffusive than the one under the white noise perturbation, which is consistent with the sharpening effect of correlated noise [50,58]. As expected, increase in noise intensity results in distributions that are more diffusive, and in distribution modes that progressively move away from the deterministic equilibrium, a phenomenon called the peak drift [50,58,59]. However, the trends in peak drift are opposite in the two perturbations; the distribution mode under white noise moves towards lower numbers of cases, while under OU noise, it shifts towards higher values.

The main difference between the two stochastic perturbations, as shown in figure 3, is the following. Whereas the distribution corresponding to OU noise stays always unimodal, the distribution corresponding to white noise undergoes a bifurcation at $\sigma \simeq 0.5\bar{\lambda }$, where an additional mode in the low COVID-19 cases regime ( $\simeq 5\times {10}^7\hspace{0.17em}\mathrm{cases}$) emerges, and the distribution becomes bimodal. This bimodality enables noise-induced transitions towards the lower number of cases, which in turn explains why the white noise model underestimates the severity of the pandemic (figure 1d).

4. Conclusion

Our study highlights the challenges of incorporating parameter uncertainties in compartmental epidemiological models. We showed that, although common, modelling these uncertainties with white noise is fundamentally flawed, since the corresponding stochastic models systematically underestimate the number of infected individuals.

We then proposed a principled modelling of the average contact rate, which accounts for the uncertainties in human social behaviour. By considering temporal correlations in each individual’s behaviour, we showed that the OU process is the correct model for uncertainties in the contact rate. We demonstrated the efficacy of our proposed stochastic model on the SIS model of sexually transmitted disease and the SEIR model of the COVID-19 pandemic.

The OU model has one unfavourable characteristic: similar to white noise, it allows for rare instances where the contact rate becomes negative, which is not epidemiologically allowed (see §S3 of the electronic supplementary material). Nonetheless, the resulting model is amenable to mathematical analysis and produces accurate forecasts that do not suffer from the idiosyncrasies of white noise. As such, the OU process serves as a reliable minimal model of the average contact rate.

Finally, the OU process has a time correlation parameter whose value needs to be derived from empirical data; such data are currently very scarce [38,60,61]. Hence, our work indicates the need for more empirical studies that quantify the correlation time in close proximity interactions.

Data accessibility

The data used in this work are available from the COVID-19 Dashboard by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University, https://github.com/CSSEGISandData/COVID-19.

The data are provided in the electronic supplementary material [62].

Authors' contributions

K.M.: data curation, formal analysis, investigation, software, visualization, writing—original draft; M.F.: conceptualization, supervision, validation, writing—review and editing.

Both authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

No funding has been received for this article.

Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.6360034.