Multi-regional COVID-19 epidemic forecast in Sweden

Abstract

The novel coronavirus disease 2019 (COVID-19) is a contagious disease with high transmissibility to spread worldwide, reported to present a certain burden on worldwide public health.

This study aimed to determine epidemic occurrence probability at any reasonable time horizon in any region of interest by applying modern novel statistical methods directly to raw clinical data. This paper describes a novel bio-system reliability approach, particularly suitable for multi-regional health and stationary environmental systems, observed over a sufficient period of time, resulting in a reliable long-term forecast of the highly pathogenic virus outbreak probability.

For this study, COVID-19 daily recorded patient numbers in most affected Sweden regions were chosen. This work aims to benchmark state-of-the-art methods, making it possible to extract necessary information from dynamically observed patient numbers while considering relevant territorial mapping.

The method proposed in this paper opens up the possibility of accurately predicting epidemic outbreak probability for multi-regional biological systems. Based on their clinical survey data, the suggested methodology can be used in various public health applications.

Key findings are:

A novel spatiotemporal health system reliability method has been developed and applied to COVID-19 epidemic data.

Accurate multi-regional epidemic occurrence prediction is made.

Epidemic threshold confidence bands given.

Introduction

The main motivation of this study was to apply a novel reliability method, newly developed by the authors, to a contemporary unfiltered health system data set.

Statistical aspects of COVID-19 (SARS-CoV-2) and other similar recent epidemics were receiving much attention in the modern research community. ¹ Generally, it is challenging to calculate realistic biological system reliability factors and outbreak probabilities under actual epidemic conditions by using conventional theoretical statistical methods.^2–11 The latter is usually due to many degrees of system freedom and random variables governing dynamic biological systems spread over extensive terrain.^12–16 In general, the reliability of a complex bio-system may be accurately assessed straightforwardly by having enough measurements or by using direct Monte Carlo simulations (if such a model is available). For COVID-19, however, the only available observation numbers are limited by the beginning of the year 2020. Motivated by the latter argument, the authors have introduced a novel reliability method for biological and health systems to predict and manage epidemic outbreaks more accurately. His study focused on COVID-19 epidemics in Sweden, focusing on cross-correlations between different regions within the same climatic zone. Sweden was chosen because of its COVID-19 origin and extensive health observations and related research available online.^4,17–29 For other studies related to statistical variations per country. ¹⁷

In this paper, an epidemic outbreak is viewed as an unexpected incident that may occur in any region of a given country at any time; therefore, the spatial spread is accounted for. Moreover, a specific non-dimensional factor $\lambda$ is introduced to predict the latter epidemic risk at any time and any place.

Biological systems are subjected to ergodic environmental influences. The other alternative is to view the process as being dependent on specific environmental parameters whose variation in time may be modelled as an ergodic process on its own.

The incidence data of COVID-19 in 21 Sweden regions from February 2020 until today were retrieved from the public website. ³⁰ As this valuable dataset is per Sweden region, the biological system under consideration can be regarded as a multi-degree of freedom (MDOF) dynamic system with highly inter-correlated regional components/dimensions. Some recent studies have already used statistical tools to predict COVID-19 development. For the linear log model, see Chu. ²

Note that while this study aims to reduce the risk of future epidemic outbreaks by predicting them, it is solely focused on daily registered patient numbers and not symptoms. For long-lasting COVID-19 symptoms, the so-called ‘long COVID’, and its risk factors and whether it is possible to predict a protracted course early in the disease. Figure 1 presents the map of Sweden's regions (counties).

Figure 1.

Left: Map of Sweden with regions (counties) with recorded COVID cases. Right: Coronavirus image.

Methods

Let one consider a MDOF bio-system represented by either response or environmental load, or combined response/load vector $\mathit{R}(t)=(X(t),Y(t),Z(t),\dots )$ , that has been either measured or simulated over a sufficiently long time period $(0,T)$ . Unidimensional $\mathit{R}$ vector component maxima being denoted as $X_T^{\max }=\underset{0\le t\le T}{max}X(t)$ , $Y_T^{\max }=\underset{0\le t\le T}{max}Y(t)$ , $Z_T^{\max }=\underset{0\le t\le T}{max}Z(t),\dots$ . Let $X_1,\dots ,X_{N_X}$ be temporally consequent local maxima of the component process $X=X(t)$ at discrete temporally increasing times $t_1^X<\dots <t_{N_X}^X$ within $(0,T)$ . Identical definitions follow for other MDOF components $Y(t),Z(t),\dots$ namely $Y_1,\dots ,Y_{N_Y};$ $Z_1,\dots ,Z_{N_Z}$ and so on. For simplicity, all $\mathit{R}(t)$ components, and hence their maxima have been assumed to be non-negative. Then

$$\begin{array}{l}P={\displaystyle \iiint {}_{\left(0,0,0,,\dots \right)}^{({\eta }_X,{\eta }_Y,{\eta }_Z,\dots )}{p}_{X_T^{\text{max}},Y_T^{\text{max}},Z_T^{\text{max}},\dots }}\\ \left(X_T^{\text{max}},Y_T^{\text{max}},Z_T^{\text{max}},\dots \right)d{X}_T^{\text{max}}d{Y}_{N_Y}^{\text{max}}d{Z}_{N_z}^{\text{max}}\dots \end{array}$$

(1)

being the probability of dynamic system survival with critical values of system components being denoted as ${\eta }_X$ , ${\eta }_Y$ , ${\eta }_Z$ ,…, $\cup$ beings logical unity operator «or» $p_{X_T^{\max },Y_T^{\max },Z_T^{\max },\dots }$ being joint probability density function (PDF) of the individual component maxima. As the system number of degrees of freedom (NDOF) is large, it is not practically feasible to estimate the joint PDF directly $p_{X_T^{\max },Y_T^{\max },Z_T^{\max },\dots }$ and therefore survival probability P. The latter probability $P$ , however, needs to be estimated, as the system expected lifetime, according to equation (1). Bio-system unidimensional components $X,Y,Z,\dots$ being now re-scaled and non-dimensionalised as follows

$$X\to {\displaystyle \frac{X}{\lambda {\eta }_X},Y\to {\displaystyle \frac{Y}{\lambda {\eta }_Y},Z\to {\displaystyle \frac{X}{\lambda {\eta }_X},\dots }}}$$

(2)

making all bio-responses non-dimensional and having the same target failure limit, when $\lambda =1$ , with target failure probability $P=P(1)$ . Equation (2) may be used now to define $P(\lambda )$ as a function of non-dimensional level $\lambda$ . Next, unidimensional system components’ local maxima being merged into one temporally non-decreasing synthetic vector $\overrightarrow{R}=(R_1,R_2,\dots ,R_N)$ in accordance with the corresponding merged time vector $t_1\le \dots \le t_N$ , $N=N_X+N_Y+N_Z+\dots$ . Each local maxima $R_j$ being actual encountered bio-system component local maxima, corresponding to either $X(t)$ or $Y(t)$ or $Z(t)$ or other system components. Constructed synthetic $\overrightarrow{R}$ vector has no data loss (see Figure 2).

Figure 2.

Example of how two components, X and Y, are merged to create a new synthetic vector $\overrightarrow{\mathit{R}}$ .

Now the non-decreasing synthetic vector $\overrightarrow{R}$ and its corresponding temporally non-decreasing occurrence times have been fully introduced.

Figure 3 presents a schematic flowchart, sketching the suggested methodology as a tool for epidemic spread surveillance.

Figure 3.

Flowchart, sketching suggested methodology.

Results

Prediction of influenza-like epidemics has long been the focus of attention in epidemiology and mathematical biology. It is well known that public health dynamics is a highly non-linear multidimensional and spatially cross-correlated dynamic system that is always challenging to analyse. Previous studies have used a variety of approaches to model influenza-like cases. This section illustrates the efficiency of the above-described methodology using the new method applied to the real-life COVID-19 datasets, presented as a new daily recorded infected patient time series spread over large terrains.

COVID-19 and influenza are contagious diseases with high transmissibility to spread worldwide with considerable morbidity and mortality. They occur most frequently seasonally in late autumn, winter and early spring, reaching their peak prevalence mostly in winter. Seasonal influenza epidemics caused by influenza A and B viruses typically occur annually during winter in temperate regions and present an enormous burden on worldwide public health, resulting in around 3–5 million cases of severe illness and 250,000–500,000 deaths worldwide each year, according to the World Health Organization (WHO). ³

This section presents a real-life application of the above-described method. The statistical data in the present section are taken from the official Sweden website. ³⁰ The website provides the number of newly diagnosed cases every day in Sweden from 22 January 2020 to 6 May 2022. Patient numbers from 21 different Sweden regions were chosen as components $X,Y,Z,\dots$ thus constituting an example of a twenty-one dimensional (21D) dynamic biological system. To unify all 21 measured time series $X,Y,Z$ ,…, the following scaling was performed according to equation (2), making all 21 responses non-dimensional and having the same failure limit equal to 1. Failure limits, or in other words, epidemic thresholds, were chosen differently for different regions in this paper ${\eta }_X,{\eta }_Y,{\eta }_Z,\dots$ and were set equal to observed 2 years maxima, twice increased. Next, all local maxima from 21 measured time series were merged into one single time series by keeping them in time non-decreasing order: $\overrightarrow{R}=(\max \{X_1,Y_1,Z_1\},\dots ,\max \{X_N,Y_N,Z_N\})$ with the whole vector $\overrightarrow{R}$ being sorted according to non-decreasing times of occurrence of these local maxima.

Figure 4 presents new daily recorded patients number plotted as a surface. Figure 5 presents the number of new daily recorded patients as a 21D vector $\overrightarrow{R}$ , consisting of assembled regional new daily patient numbers. Note that vector $\overrightarrow{R}$ does not have physical meaning on its own, as it is assembled of different regional components with different epidemic backgrounds. Index j is just a running index of local maxima encountered in a non-decreasing time sequence.

Figure 4.

New daily recorded patients number plotted as a surface: Provinces and time dependency.

Figure 5.

Number of new daily recorded patients as 21D vector $\overrightarrow{\mathit{R}}$ . Left: As it is. Right: Scaled by equation (12).

Figure 6 (left) presents 100 years return level extrapolation according to equation (9) towards epidemic outbreak with 100 year return period, indicated by the horizontal dotted line, and somewhat beyond, $\lambda =0.1$ cut-on value was used. Dotted lines indicate extrapolated 95% confidence interval according to equation (10). According to equation (5), $p(\lambda )$ is directly related to the target failure probability $1-P$ from equation (1). Therefore, in agreement with equation (5), system failure probability $1-P\approx 1-P_k(1)$ can be estimated. Note that in equation (5), N corresponds to the total number of local maxima in the unified response vector $\overrightarrow{R}$ . Conditioning parameter $k=5$ was found to be sufficient due to occurrence of convergence with respect to k, see equation (6). Figure 6 exhibits reasonably narrow 95% CI. The latter is an advantage of the proposed method.^31–35

Figure 6.

Left: 100 years return level (horizontal dotted line) extrapolation of ${\mathit{p}}_{\mathit{k}}(\mathit{\lambda })$ towards critical level (indicated by a star) and beyond. Extrapolated 95% CI indicated by dotted lines. Right: COVID Swedish national statistics, SODP plot.

Note that while being novel, the above-described methodology has a clear advantage of utilising available measured data sets quite efficiently due to its ability to treat health system multi-dimensionality and perform accurate extrapolation based on a quite limited data set. Note that the predicted non-dimensional $\lambda$ level, indicated by a star in Figure 6 (left), represents the probability of an epidemic outbreak in any Sweden region in the years to come.

The second-order difference plot (SODP) originated from the Poincare plot. SODP provides observing the statistical situation of consecutive differences in time series data.

Figure 6 (right) presents the SODP plot. This kind of plot can be used for data pattern recognition and comparison with other data sets, for example, for the entropy artificial intelligence (AI) recognition approach. ³⁶ Note that EVT is asymptotic and 1DOF, while this study introduces MDOF and sub-asymptotic approaches. To summarise, the predicted non-dimensional λ level, indicated by the star in Figure 6 (left) represents the probability of world cancer deaths in the years to come. The methodology's limitation lies in its assumption of the underlying bio-environmental process quasi-stationarity.^37,38

Discussion

Traditional health systems reliability methods dealing with observed time series do not have the advantage of dealing efficiently with systems possessing high dimensionality and cross-correlation between different system responses. The key advantage of the introduced methodology is its ability to study the reliability of high dimensional non-linear dynamic systems.

Despite the simplicity, the present study successfully offers a novel multidimensional modelling strategy and a methodological avenue to implement the forecasting of an epidemic during its course.

This paper studied recorded COVID-19 patient numbers from 21 different Sweden regions, constituting an example of a 21 dimensional (21D) observed in 2020−2022. The novel reliability method was applied to new daily patient numbers as a multidimensional system in real-time. The theoretical reasoning behind the proposed method is given in detail. Note that the use of direct either measurement or Monte Carlo simulation for dynamic biological system reliability analysis is attractive; however, dynamic system complexity and its high dimensionality require the development of novel robust and accurate techniques that can deal with a limited data set at hand, utilising available data as efficient as possible.

The main conclusion is that Sweden's public health system under local environmental and epidemiologic conditions is well managed. The predicted 100-year return period risk level $\lambda$ of an epidemic outbreak is very low.

Various authors with different approaches have shown the usage of statistics through EVT and other models in medicine. One such method used the block maxima approach, while another used the Peak Over Threshold approach to estimate the distribution of extremes. Even though both these studies showed their suitability for estimating the extreme values, each of them had its limitations, with one of them requiring a large amount of data.

This study further aimed to develop a general-purpose, robust, and straightforward multidimensional reliability method. The method introduced in this paper has been previously validated by application to a wide range of simulation models, but for only one-dimensional system responses and, in general, very accurate predictions were obtained. Both measured and numerically simulated time series responses can be analysed. It is shown that the proposed method produced a reasonable confidence interval. Thus, the suggested methodology may become appropriate for various non-linear dynamic biological systems reliability studies. Finally, the suggested methodology can be used in many public health applications. The presented COVID-19 example does not limit areas of new method applicability.

The major limitations of the suggested approach are as follows:

Manipulated underlying data set

Underlying trend

System non-stationarity

While, as mentioned, C is not a major obstacle for the suggested method, future works may extend current findings to cope with B's limitation, namely to identify the underlying trend.