A stochastic SEIHR model for COVID-19 data fluctuations.

Although deterministic compartmental models are useful for predicting the general trend of a disease's spread, they are unable to describe the random daily fluctuations in the number of new infections and hospitalizations, which is crucial in determining the necessary healthcare capacity for a specified level of risk. In this paper, we propose a stochastic SEIHR (sSEIHR) model to describe such random fluctuations and provide sufficient conditions for stochastic stability of the disease-free equilibrium, based on the basic reproduction number that we estimated. Our extensive numerical results demonstrate strong threshold behavior near the estimated basic reproduction number, suggesting that the necessary conditions for stochastic stability are close to the sufficient conditions derived. Furthermore, we found that increasing the noise level slightly reduces the final proportion of infected individuals. In addition, we analyze COVID-19 data from various regions worldwide and demonstrate that by changing only a few parameter values, our sSEIHR model can accurately describe both the general trend and the random fluctuations in the number of daily new cases in each region, allowing governments and hospitals to make more accurate caseload predictions using fewer compartments and parameters than other comparable stochastic compartmental models.

Introduction

Over the past year, the coronavirus disease 2019 (COVID-19) pandemic has placed enormous stress on healthcare systems worldwide. The severity of the pandemic, given the fact that it takes a long time to develop effective vaccines, has imposed tremendous pressures and responsibilities onto the healthcare systems in all countries and regions, especially those with limited medical resources such as available staff, equipment and facilities. In particular, differing from most infectious diseases found to date, COVID-19 is especially violent, aggressive and fast-spreading, and even among those deemed “recovered” from the disease, there are many for which adverse effects have lingered for months after the initial symptoms [1]. Meanwhile, failure to meet the demand for hospital resources can lead to resource saturation and a growing backlog of infectious patients requiring hospitalization, in turn increasing the total transmission rate due to unisolated infectious individuals and causing an adverse feedback loop.

Related studies on COVID-19

Due to an incomplete knowledge of COVID-19 during the early stages, scientists used existing models to forecast the pandemic and made inaccurate predictions. Fortunately, based on historical data, people can still make relatively accurate predictions by applying some kinds of model-free methods [2, 3]. However, to generate policy-relevant insights into the nonpharmaceutical interventions, people still need to understand the physical principles of the pandemic. Since the clinical data of COVID-19 have been carefully studied  [4], scientists now can make more accurate epidemic predictions using corresponding mathematical models [5, 6]. Reference [7] presents the details of three regional-scale models for predicting and estimating the trend of the pandemic. We also introduced a compartmental model and made an accurate trend forecasting on Hong Kong’s COVID-19 pandemic  [8].

To predict the number of cases in a region and the corresponding healthcare demands, epidemiological models are very useful and have been constructed to describe the spread of diseases like COVID-19. In particular, compartmental models, which track the evolution of the number or proportion of the individuals in each disease state (susceptible, infected, removed, etc.) using a system of ordinary or stochastic differential equations, have been developed for nearly a century and have proved effective and useful [9, 10].

Contributions of this paper

In this paper, we develop a stochastic version of the modified SEIHR compartmental model established in [8], which allows one to accurately model the random fluctuations in the number of daily new cases in the COVID-19 epidemic process. We prove sufficient conditions on the parameters of the new stochastic SEIHR (sSEIHR) model, such that the disease-free equilibrium (DFE) of the system is stochastically stable, based on an estimated basic reproduction number. Our extensive numerical results demonstrate strong threshold behavior near the estimated basic reproduction number, suggesting that the necessary conditions for stochastic stability are actually close to the sufficient conditions established.

The proposed sSEIHR model is then fitted to public COVID-19 data from various regions worldwide, and it is demonstrated numerically that the model can accurately describe both the general trend and the random fluctuation in the number of daily new cases in each region. The ability of our model to accommodate a wide range of geographical regions and the variations in their various COVID-19 outcomes with minimal adjustments to model parameters show the robustness of our model. This can offer suggestions to governments and hospitals for their making more accurate caseload predictions compared to using a deterministic model.

Recall that it has been demonstrated in [8] that deterministic models using fewer compartments than the SEIHR model, e.g., SIR or SEIR, fail to accurately describe the evolution of COVID-19. Therefore, except potentially by means of more complex compartment interactions, both the SEIHR model developed in [8] and the new sSEIHR model developed here are the smallest compartmental models that can accurately describe the number of active COVID-19 cases in a region of concern. The model’s performance suggests that it captures key hidden dynamical features of general epidemics in general environments and deserves independent investigation.

The true significance of this work is the generality of the proposed model with few parameters. In particular, it is surprising to find that the model fits very well to a number of regions/countries worldwide with only a few parameter value changes, under the consideration of the heterogeneity of economic and social features, medical resources, quarantine measures, mitigation strategies and healthcare systems’ characteristics in different regions/countries. In addition, within each region, changes in non-pharmaceutical intervention (NPI) strength over time can be expressed as a single parameter, while all other parameters remaining constant. The robustness and accuracy observed make this model valuable for real applications.

Background and related work

Notation

For vectors , let denote the set and denote the set .

The SIR model

A classical approach for epidemiological modeling is to use compartmental models [9]. One such compartmental model is the SIR model, in which individuals are classified as susceptible (S), infected (I), or removed (R) (either recovered or deceased). The number of individuals in each state, i.e., compartment, can be described by a system of ordinary differential equations.

Let , , and denote the numbers of susceptible, infected, and removed individuals, respectively, and let N be the population size, at time . The system can then be described bywithwhere , with and . For this model, the basic reproduction number is important, which can be thought of as the number of expected cases directly generated by an infected individual, if all other individuals are susceptible to infection. The epidemic is expected to persist if , but will die out if .

It was proved [11] that for , uniquely determines whether the epidemic through system (1) will die out () or persist (), independent of the nonzero initial conditions. This threshold behavior demonstrates the power and usefulness of compartmental modeling in the study of epidemic spreading.

Other deterministic models

Numerical results shown in [8] demonstrate that the SIR model is inadequate for describing COVID-19 dynamics. Reasons include:To account for the above issues simultaneously, in [14] a nine-compartment model is constructed, and the model is referred to as -SEIHRD, in which there are two I, H, and R compartments, respectively. However, while the -SEIHRD model is useful for detailed study of COVID-19 dynamics, it may not be desirable to use such a complicated model for more effective parameter fitting, especially if only the total number of active cases is required to be estimated. This problem is resolved in [8], with a model having only five compartments:with , , and , where S, E, I, H, and R denote the proportions of susceptible, exposed, infectious, hospitalized, and removed individuals, respectively. Note that in our model, “infectious” denotes symptomatic infectious individuals, whereas “exposed” includes both non-infectious exposed individuals and asymptomatic infectious individuals. The parameter allows for asymptomatic transmissions, a key feature of COVID-19 propagation. The meaning of each parameter in (2) is given in Table 1. Further interpretation about the compartments can be found in  [8].

Table 1

Meaning of parameters in the SEIHR model (2)

Parameter	Definition
η	Transmission rate of exposed individuals
α	Transmission rate of (symptomatic) infected individuals
β	Reciprocal of the mean latent period, i.e., the rate at which exposed individuals become symptomatic
γ	Rate at which infected individuals are hospitalized
ω_E	Rate of recovery of non-hospitalized exposed individuals
ω_I	Rate of recovery of non-hospitalized infected individuals
ω_H	Rate of recovery of hospitalized individuals

Like some other diseases, COVID-19 has a long latent period. The SIR model does not take this into account, however. As an improvement, SEIR models have been developed to model the effect of latent periods, where E represents an additional “exposed” compartment [12].

Ideally, symptomatic individuals are detected and treated in hospitals, wherein they cannot infect the general community. The SIR model does not take this into account, either. As an improvement, the SEIQR or SEIHR model creates an additional compartment for such individuals, with Q and H representing “quarantined” and “hospitalized”, respectively [13].

The COVID-19 pandemic has a large number of asymptomatic infections. The SIR model certainly does not take this into account. Therefore, some other models seek to divide infected individuals into symptomatic and asymptomatic compartments. Alternatively, the E compartment in the SEIR model can be made infectious [8].

The basic reproduction number of system (2) isHereafter, the term “deterministic SEIHR model” or “dSEIHR model” refers to the one described by system (2). A diagram showing the transitions of the dSEIHR model (2) is given in Fig. 1. Note that (2) differs from the traditional SEIHR model in that the exposed individuals may also be infectious.

Fig. 1

Graphical depiction of the dSEIHR model

To further improve the model ability in prediction, especially over a long time period, a factor was introduced in [8] to denote the strength of various interventions (a smaller p(t) implies a stronger intervention), leading towhich is called the p-dSEIHR model.

For fitting the p-dSEIHR model (3) to real COVID-19 data, piecewise constant functions were used in [8] to approximate . This matches the discrete nature of governmental and healthcare policy, where occasional policy changes will cause immediate changes in .

Stochastic epidemic models

While the deterministic models described in Sects. 2.1 and 2.2 have proven useful for estimating the cumulative number of infections over time, these do not capture the stochastic nature of disease propagation. This can lead to healthcare dimensioning and resource allocation problems if the actual number of daily new cases is significantly different from the predicted value, which in turn can lead to resource saturation and a growing backlog of infectious patients requiring hospitalization, increasing the total transmission rate due to unisolated infectious individuals and causing an adverse feedback loop.

To describe random fluctuations, stochastic compartmental models are commonly used, where the system in Sect. 2.2 is modified to be a drift-diffusion process of the form:where f and g are vector-valued functions and denotes a Wiener process, i.e., derived by integrating Gaussian white noise.

Examples of stochastic compartmental models in the literature include [15] and [16], which are the stochastic analogs of the deterministic SIR and -SEIHRD models described in Sect. 2.2, respectively. Notably, in [15, 16] it is assumed that the noise in these two stochastic models is attributable to a single state transition, i.e., in [15] and in [16]. The same assumption will be adopted for the new sSEIHR model.

Stochastic SEIHR model with Gaussian white noise

Model formulation

While the deterministic SEIHR model (2) addresses the COVID-19 phenomena listed in Sect. 2.2, it cannot model the random fluctuations in the number of daily new cases. To address this issue, we therefore convert system (2) to a stochastic SEIHR (sSEIHR) model by introducing an additive random component to the transition, thus obtaining the following drift-diffusion process:with and , where W(t) denotes a Wiener process obtained from integrating Gaussian-distributed white noise. Here, in this model, only the two terms I and H are considered having random fluctuations due, for example, to uncertain hospital management.

A graphical depiction of system (5) is shown in Fig. 2, where denotes Gaussian white noise.

Fig. 2

Graphical depiction of the sSEIHR model (5), where denotes Gaussian white noise

Stability results

It is easy to see that (5) has a disease-free equilibrium . Regarding the system stability at this equilibrium, we have the following main result.

Theorem 1

Ifthen the DFE of system (5) is stochastically stable, in the sense defined in “Appendix A.”

For a proof of Theorem 1, see “Appendix B”.

Note that in the absence of noise, i.e., , becomes the basic reproduction number of the corresponding deterministic system (2).

Numerical results

Consider a system with a total population of , where , , , , and denote the numbers of susceptible, exposed, infected, hospitalized and removed individuals, respectively, with , , , , and , which evolve according to the stochastic differential Eq. (5).

The impact of the noise term on the epidemic process is examined, using the Euler–Maruyama method with initial state , , and , with the system parameters , , and .

Figure 3 shows the simulated value of and , for various values of and . It shows that the value of undergoes a sharp transition in the vicinity of the line , suggesting that the necessary conditions for stochastic stability are close to the sufficient conditions provided in Theorem 1. Furthermore, Fig. 3 demonstrates that when the noise parameter increases, the final value slightly decreases.

Fig. 3

Comparison between the simulated value of and of the sSEIHR model (5) for the initial values and parameter settings defined in Sect. 3.3

Figure 4 shows the impact of the noise parameter on the epidemic process. Note that for large values of , the probability becomes significant, such that patients can return to the state I from state H and even cause H(t) to become negative, as can be seen clearly from Fig. 4a. Therefore, smaller values of are generally more suitable for modeling the random fluctuations at the hospitalization rate of the system. Furthermore, as the primarily concern is to estimate some rough probabilistic bounds on the number of daily new cases, the simulation with Gaussian white noise is sufficient.

Fig. 4

Impact of the noise intensity on the epidemic process described by the sSEIHR model (5)

Real-data analysis

Analogous to our modification of the dSEIHR model to produce the p-dSEIHR model, the sSEIHR model is modified by introducing a scaling factor p(t) to denote the strength of various anti-epidemic interventions:This will be referred to as the p-sSEIHR model.

Hong Kong

Now, consider the number of daily new cases of COVID–19 in Hong Kong, between June 25 and October 8, 2019. This is compared to the expected values obtained from the p-dSEIHR model (3).

The fitted parameters of the p-dSEIHR model are given in Table 2, where and for t in each period shown, with for t prior to . Note that this dataset corresponds to and H is the only observable compartment in the model. Note that the parameters N, , , , and () are fixed parameters, so only , , and () are fitted.

Table 2

Parameter setting for the Hong Kong dataset

Parameter	Value	Parameter	Value
N	7 500 700	p_1	0.635
φ	0.03	t_1	4/7∼23/7
η	0.48	p_2	0.345
α	0.5	t_2	24/7∼13/8
β	0.14	p_3	0.235
γ	0.7	t_3	14/8∼27/8
ω	0.1	p_4	0.4445
		t_4	28/8∼8/10

Since the noise term in the p-sSEIHR model is proportional to I(t), we are interested in the values . As shown in Fig. 5b, the values roughly follow a left-truncated Gaussian distribution with mean and standard deviation . Figure 5a shows , , and for the Hong Kong dataset.

Fig. 5

Ratios and the estimated -range of fluctuations in the number of daily new COVID-19 cases in Hong Kong, where refers to 25 June 2020

Next, using the same fitted parameters as that in the p-dSEIHR model, the task is to find a value for such that has an standard deviation of , where denotes the mean value of obtained from 1000 sample paths of the p-sSEIHR model. From the data, the result is found to be .

Figure 6 shows the cumulative and daily numbers of COVID-19 cases in Hong Kong between June 25 and October 8, 2019, the p-dSEIHR model estimate, and the (2.5,97.5)-percentile range of sample paths from the corresponding p-sSEIHR model. As shown in Fig. 6b, this range (colored gray) closely matches the envelope obtained from the p-dSEIHR model for all time points.

Fig. 6

Actual and predicted numbers of COVID-19 cases in Hong Kong, as estimated using the p-sSEIHR model (5)

Other regions worldwide

In simulation, the process described in Sect. 3.4.1 is repeated for four global regions, namely Germany, Spain, South Africa, and the state of New York, using data obtained from [17]. The parameter settings of the fitted p-dSEIHR and p-sSEIHR models are given in Table 3. Again, the values of roughly follow a left-truncated Gaussian distribution, with standard deviations ranging from 0.24 to 0.33.

Table 3

Parameter settings for four global regions, with , , and for all regions

	Germany	Spain	South Africa	New York
N	83 783 945	47 431 256	59 622 350	19 745 289
φ	0.035	0.035	0.03	0.03
η	0.338	0.57	0.47	0.638
α	0.46	0.6	0.54	0.7
p_1	0.476	0.635	0.47	0.87
t_1	27/3∼25/5	13/3∼27/3	22/3∼16/4	13/3∼21/3
p_2	0.64	0.345	0.555	0.35
t_2	26/5∼8/7	28/3∼12/4	17/4∼12/7	22/3∼14/4
p_3	–	0.235	0.378	0.253
t_3	–	13/4∼10/5	13/7∼13/9	15/4∼14/5
p_4	–	0.4445	0.47	0.31
t_4	–	11/5∼31/5	14/9∼12/10	15/5∼6/7

Figure 7 shows the distributions of for each of the four global regions, as well and the mean and standard deviation for each distribution.

Fig. 7

Distributions of in four global regions

Figure 8 shows the cumulative and daily numbers of COVID-19 cases in each global region during the specified time periods, the corresponding p-dSEIHR estimates, and the (2.5,97.5)-percentile ranges of sample paths from the corresponding p-sSEIHR model. It can be seen that these ranges, colored gray, closely match the envelope obtained from the p-dSEIHR model for the number of daily new cases (right column) in each region. This result not only demonstrates the usefulness of the p-sSEIHR model, but also reveals that the noise profiles in different regions may have a similar nature.

Fig. 8

Real COVID-19 data and p-sSEIHR predicted values for four global regions

Discussion

Section 3.2 shows the theoretical significance of the sSEIHR model. By studying the stability of the disease-free equilibrium, an estimated basic reproduction number was obtained, which helped understand the critical condition of the COVID-19 outbreak. Section 3.3 confirms the numerical results in accordance with the theoretical results. The numerical results also show the critical value of the noise intensity where the pandemic happened. Most importantly, Sect. 3.4 shows the accuracy of our sSEIHR model on COVID-19 pandemic predictions in different irrelevant regions.

The key to the success achieved on COVID-19 pandemic forecasting/fitting can be appreciated by as follows. First, our deterministic SEIHR model captures the critical features of the COVID-19 transmission process, such as the disease transmission during the latent period, the hidden infections, and the isolation policy. Then, using only a few parameters, we successfully developed a deterministic model that accurately described the COVID-19 pandemic trend. In order to predict the daily caseload, the corresponding stochastic SEIHR model only needs one more parameter than the deterministic SEIHR model. Surprisingly, the real data analysis shows that the noise intensities in different regions are very close to each other. In the end, we used one parameter to describe the policy changes during the pandemic. Based on the timing of local intervention in each region, we changed the parameter accordingly and obtain an accurate prediction of the trend of the COVID-19 pandemic. All of this proves that our model has excellent performance in COVID-19 pandemic prediction.

Concluding remarks

In this paper, we developed a model described by a system of stochastic differential equations to describe a stochastic SEIHR (sSEIHR) process. We proved that the disease-free equilibrium is stochastically stable when the estimated basic reproduction number is less than one. Extensive numerical results demonstrate strong threshold behavior near the estimated basic reproduction number, suggesting that the necessary conditions for stochastic stability are close to the sufficient conditions given in Theorem 1. We also included a factor p(t), denoting the strength of various anti-epidemic interventions, to establish a p-sSEIHR model, and used this model to fit the daily number of hospitalized cases in Hong Kong and other four global regions worldwide. We find that the five regions have a very similar noise profile, with an estimated value from 0.24 to 0.33 and a fitted parameter from 0.03 to 0.035 for the p-sSEIHR model.

As the new p-sSEIHR model not only predicts the mean number of daily new cases of COVID-19 in a region, but also specifies a possible range of values, it can provide healthcare practitioners with a clearer picture of potential demands for medical resources due to epidemics. Furthermore, the similarity of the noise parameter for different regions suggests that the information about a disease learned from one region can be used for other regions regarding management policies, even at a less advanced stage of an epidemic, as long as the transmission parameters and in the model can be estimated with a reasonable accuracy.