Mathematical analysis for COVID-19 resurgence in the contaminated environment

1.   Introduction

COVID-19 is an epidemic of novel coronavirus pneumonia, which is enveloped, single-stranded, positive-sense RNA virus belonging to the family of Coronaviridae ^[1]. The epidemic, nowadays, has spread all over the world. The clinical symptoms are fever, dry cough and fatigue, and patients may develop into severe states ^[2]. The epidemiological investigations indicate that the epidemic disease results from direct transmissions by contact between susceptible individual and infectious patient, or from indirect transmission where susceptible individuals are exposed to pathogens in the contaminated environment ^[3]. Close contact with infectious individuals is the main channel of direct transmission. According to WHO's guidance document for COVID-19 ^[4], infected individuals can carry or produce virus. Susceptible population may inhale the respiratory droplets, shedding from infected individuals, and then may become exposed or asymptomatic individuals, which depends on personal situation ^[3]. Asymptomatic individuals constitute a special group of infectious patients who have no clinical symptoms ^[3,5], and their transmissions are harder to be noticed. Both of exposed individual and asymptomatic individuals have the ability to transmit COVID-19 by talking together, dining together and living together ^[5,6].

Though the direct transmission is the main route of COVID-19 transmission, many infection cases caused by indirect transmission have already occurred in China. The main channels of indirect transmission are inhaling contaminated aerosols and contacting public facilities with the virus. The primary infection between the individuals living on adjacent floors in Haerbin of China in April was probably infected by touching buttons with the virus or inhaling contaminated aerosols in the elevator ^[7]. In addition, the infection of a laundrywoman occurring in Jilin of China in May had a great probability of indirect transmission by touching the clothes with the virus ^[8]. Moreover, the local prevalence of COVID-19 in Beijing of China in June originated from indirect transmission by touching the contaminated frozen meat in the Xinfadi wholesale market. Since the coronavirus can survive in the environment for several days ^[9], it is important to investigate how the indirect transmission contributes to the outbreak of COVID-19.

A number of dynamical models for COVID-19 are proposed to depict the patterns of disease progression. By comparing with SARS and MERS, Jin et al. investigate the critical parameters of the epidemic in Hubei of China, estimate the efficacy of control measures and give helpful suggestions ^[10]. Tang et al. propose a compartmental model, estimate the control reproduction number, and assess the cost and efficacy of epidemic prevention to provide the control policy for reopening time ^[11,12]. From the early data in China, Zhao et al. estimate the basic reproduction number of China ^[13], and Hu et al. estimate the basic reproduction numbers of COVID-19 in different regions of China and emphasize the transmission risk of input infections ^[14].

Notice that most studies by compartmental modeling focus upon the direct transmissions from close contacts. The objective of this paper is to propose a mathematical model of COVID-19 that includes not only the direct transmission by close contacts, but also the indirect transmission from environmental pathogens. Because patients with clinical symptoms are quickly quarantined to the hospital, we will focus our attention on influences of exposed and asymptomatic infections on the disease dynamics. By mathematical analysis, we find the relationships among the basic reproduction number, peak value and final size. By numerical simulations, we find that the input of asymptomatic individuals is dangerous to the epidemic spread. Furthermore, the detection intensity and environmental disinfection significantly affect the basic reproduction number and total infection size. Hence, the epidemic peak, epidemic time and final infection size can be effectively limited through these prevention measures.

The organization of remaining paper is set as follows. The model formulation and analysis are given in Section 2. Numerical simulations are implemented in Section 3 to extrapolate how indirect transmissions and covert infections contribute to the resurgence of COVID-19 in Heilongjiang of China. In the end, we present discussions and draw some conclusions.

2.   Model formulation and analysis

2.1.   Model formulation

We start from the formulation of a mathematical model. The population is divided into five compartments: susceptible, asymptomatic, exposed, hospitalized and recovery individuals, the numbers of which at time $t$ are denoted by $S(t), A(t), E(t), H(t)$ and $R(t)$ respectively. The compartment of infectious individuals with clinical symptoms is neglected since the detection measures are implemented so that they are quickly quarantined to the hospital. We differentiate the asymptomatic individuals from exposed individuals because there are differences between asymptomatic individuals and exposed individuals. Indeed, an asymptomatic individual has a probability of self-healing, and has a longer duration than an exposed individual ^[15].

Note that an exposed individual has the potential to infect others, and so does an asymptomatic individual, which has been confirmed in epidemiological investigations ^[3,6,7,9]. For example, the infection cases caused by exposed individuals account for $6.30\%$ in their close contacts, which is slightly bigger than $4.11\%$ for asymptomatic transmissions in Ningbo of China ^[6]. Let $\beta_1$ and $\beta_2$ denote the valid contact coefficients of exposed and asymptomatic individuals through direct transmission respectively. Then the incidence from direct transmissions is described by

We now consider the indirect transmissions which are from aerosols or the public facilities with the virus. According to medical investigations, coronavirus can survive for several days in damp, dark and cold environment ^[7,8,9]. For the people who stay in such a contaminated environment for a longer time $\tau$ because of service duty or other reasons, the accumulated effect must be considered. Let $V(t)$ be the concentration of virus in environment at time $t$ and $\beta_3$ be the rate at which a susceptible individual takes in virus from the contaminated environment. This could be the case by touching the pathogens through hands or clothes while he or she stays in the contaminated environment. If $\eta$ is the decay rate of virus with a host, then the accumulated virus concentration at time $t$ within the susceptible individual is

As a result, one can assume that the infection rate due to indirect transmission at time $t$ is described by

For convenience in notation, we denote the total infection force of disease transmission by

The susceptible individuals after infected, go either into the asymptomatic compartment at proportion $p$, or into the exposed compartment at proportion $1-p$, where $p\in[0, 1]$. Assume that the exposed and asymptomatic individuals transit into the hospital at rates $\sigma$ and $\alpha$ respectively. Moreover, asymptomatic members and hospitalized individuals are recovered at rate $\gamma_1$ and rate $\gamma_2$ respectively. Furthermore, we let $\xi_1$ and $\xi_2$ denote the production rates of virus from asymptomatic and exposed individuals, respectively. In addition, since the average death rate of COVID-19 patients in China is 2.1% ^[16] and 1.37% in Heilongjiang province, China ^[7], which is quite low, we neglect the impact of fatality. The schematic diagram of disease transmission and progression is shown in Figure 1, and the dynamics of the state variables are described by the following differential equations:

where a prime denotes the differentiation of a state variable with respect to time $t$. Since the first four equations in system (2.1) are independent of $H$ and $R$, it suffices to consider following system

with initial values

Since there is no available data to show the distinction between $\beta_1$ and $\beta_2$ and the distinction between $\xi_1$ and $\xi_2$ in current medical investigations, we confine ourselves to $\beta_1 = \beta_2 = \beta$ and $\xi_1 = \xi_2 = \xi$ in the present paper.

2.2.   The basic reproduction number

The basic reproduction number is an average number of secondary infections produced by a single infection in an entirely susceptible population during its infection period. By the next generation matrix method ^{[17,18,19,20]}, the last three equations of system (2.2) can be written as

where $x = (A, E, V)^T$, and

Set

Through the method in Zhao ^[17], the basic reproduction number $\mathcal{R}_0$ of system (2.2) is given by $\rho(FU^{-1})$, the spectral radius of $FU^{-1}$. Thus,

where

2.3.   Analysis of peak value and final size

In order to qualitatively investigate the peak value of infections, following Feng ^[21], we introduce a weighted sum of infectious individuals and pathogens:

Clearly, $Y(t)$ is an index of disease infection at time $t$ because the infection risk becomes higher with the increase of $Y$.

Define the control reproduction number by

We have the following relation between $\mathcal{R}_c$ and $\mathcal{R}_0$.

Theorem 2.1. ${\rm{sign}} (\mathcal{R}_0-1) = {\rm{sign}} (\mathcal{R}_c-1)$.

Proof. Note that $\mathcal{R}_0$ and $\mathcal{R}_c$ can be written as

where

Theorem 2.1 is evident if $a\geq 1$. As $a < 1$, it is easy to see that

is equivalent to

Therefore,

The following theorem indicates that the existence of a peak for the infection index $Y$ is determined by the control reproduction number $\mathcal{R}_c$.

Theorem 2.2. The infection index $Y$ hits a peak $Y_{\rm max} = Y(t^*)$ for some $t^* > 0$ if $\mathcal{R}_c > 1$, and has no peak if $\mathcal{R}_c < 1.$

Proof. First, we claim $Y(t)\to 0$ as $t\to \infty$. To this end, we let

where $a > \max\{\xi/(\alpha+\gamma_1), \xi/\sigma\}$. Then we get $U'(t) \leq 0$. In the set $\{(S, A, E, V):U' = 0\}$, it is easy to see $S$ axis is the largest invariant set. From the LaSalle's invariance principle, we conclude

This proves the claim.

Based on the second, third and fourth equations in system (2.2), we get

where $\theta \in (0, \tau)$ and

If $\mathcal{R}_c > 1$, we can easily get $Y'(0) > 0$. Since $\lim\limits_{t\to+\infty}{Y(t)} = 0$, there must be a $t_1 > 0$ such that $Y_{\rm max} = Y(t_1)$. If $\mathcal{R}_c < 1$, we have $Y'(0) < 0$. Since the subsystem that consists of the last three equation of system (2.2) with a decreasing function $S(t)$ is a monotone system, we see that $Y'(t) < 0$ for $t > 0$. Therefore, $Y(t)$ has no peak for $t > 0$.

Final size $z$ is a quantitative value which embodies the influence of epidemic ^[22,23,24], and is defined as

We state that the final size $z$ of system (2.2) satisfies

where

Theorem 2.3. The final size $z$ of system (2.2) is the unique root of Eq (2.4).

Proof. By combining $S$ with $A$ and $E$, we get following two equations

From the third equation of system (2.2) and Eq (2.5), we get the following three equations

where

$(S_0, A_0, E_0, V_0)$ and $(A_{\infty}, S_{\infty}, E_{\infty}, V_{\infty})$ are initial and terminal values of $(S, A, E, V)$, respectively. Solving Eq (2.6) and using Eq (2.3), we obtain

Next, through integral transformation $t-r = s$, we get the following integration

By dividing $S$ in the first equation of system (2.1) and then integrating, we have

Using Eqs (2.7), (2.8) and (2.9), we yield

Further, we have

Finally, we get

where

Since $f(0)*f(1) < 0$, it is easy to see that (2.10) admits a root $z_0$. To examine the uniqueness, we let

Then $f'(z_{\min})\leq0$. If there exists another root $\bar z$, it follows from Rolle's theorem that there exists a $z_\xi\in(z_{\min}, \bar{z})$ such that

But $f'(z) = -1+\mathcal{R}_c \exp\{-\mathcal{R}_cz-\zeta\}$ is a strictly monotone decreasing function. Therefore, Eq (2.11) cannot hold. Consequently, we conclude that Eq (2.10) admits a unique root.

The final size $z$ can only be numerically solved. However, it paves a way to estimate the influences of parameters on the disease spread. In particular, one can infer how different transmission routes contribute to the disease outbreak. For the latter purpose, we take $p = 0$ in Eq (2.4) to obtain the formula for the infection size from exposed individuals, and $p = 1$ for the infection size from asymptomatic individuals. We also set $\xi = 0$ to see how indirect transmissions contribute to the final size of epidemic spread.

As $p = 0$, from Eq (2.4) we get the final size equation

where

and $\mathcal{R}_c$ is

As $p = 1$, Eq (2.4) implies that the final size equation is given by

where

and $\mathcal{R}_c$ satisfies

As $\xi = 0$, from Eq (2.4), we obtain the final size equation

where

and $\mathcal{R}_c$ is described by

In summary, through the mathematical analysis in this subsection, we have obtained the computation formulae for the basic reproduction number and control reproduction number, and established the equations to determine the final sizes of COVID-19 spread in several important scenarios. These provide the basis to calculate the characteristics of COVID-19 resurgence from real data in the next section.

3.   Numerical simulation

3.1.   Raw data

In April of 2020, COVID-19 epidemic resurged in Heilongjiang province of China. Health Commission of Heilongjiang Province provided the number of asymptomatic infections per day, where asymptomatic individuals are those whose nucleic acid tests are positive but do not show any symptoms of infection ^[16]. We acquire the real data from Heilongjiang Provincial Health Committee ^[7]. Table 1 shows the number of daily asymptomatic individuals and Table 2 gives the number of cumulative confirmed cases. The histograms of asymptomatic individuals and cumulative confirmed individuals of Heilongjiang are shown in Figure 2. The number of asymptomatic individuals in Heilongjiang kept increasing in the first week in April and hit a peak on April 7th. It began to decline slowly from April 8th and got extinct at the end of April with the final size of 460 confirmed cases.

3.2.   Parameter estimation

Previous studies in the transmission of COVID-19 have given the information for most of parameter values in system (2.2). Specifically, the natural latent period of asymptomatic individuals is estimated between $3-14\ {\rm days}$ ^[15,25], and the natural latent period of exposed individuals is estimated between $1-6\; {\rm days}$ ^[11,15]. The recovery rates of asymptomatic individuals is estimated at $1/12\; {\rm day}^{-1}$, and the recovery rate of hospitalized individuals is estimated at $1/10\; {\rm day}^{-1}$ ^[16,26]. Environmental coronavirus has its survival duration that ranges from $3$ hours to $3\; {\rm days}$ ^[9]. The probability of being asymptomatic is estimated between $10\%-30\%$ or even higher in some clustering infection cases ^[5,27]. Moreover, we get from Heilongjiang province people's government ^[28] that its total population is around $3.7*10^7$ individuals, which is fixed as the initial size of susceptible. Further, we obtain from ^[7] that the initial size of exposed individuals is 22, and the initial size of asymptomatic individuals ranges from 15 to 39.

What we need to estimate are the valid transmission coefficients $\beta, \beta_3$, the virus production rate $\xi$, the virus decay rate $\eta$, and the initial concentration $V_0$ of environmental virus. Through the least square method for the data from Table 1 and Table 2, we obtain $\beta = 7.8*10^{-9}\ \text{person}^{-1}\text{day}^{-1},$ $\beta_3 = 9.56*10^{-14}\ \text{person}^{-1}\text{day}^{-1}$, $\eta = 0.63\; {\rm day}^{-1},$ $\xi = 1.12*10^2\; \text{person}^{-1}\text{day}^{-1}\text{ L}^{-1}$ and $V_0 = 9.65*10^6/{\rm L}$. Consequently, we can present Table 3 for the estimated values or ranges of the parameters and initial values in system (2.2). Next, we choose

which fall in their ranges as indicated in Table 3. We show in Figure 3 that the predicated values by model (2.2) fit the real data very well. Further calculations show that the peak values $(A_{\rm peak}, E_{\rm peak}) = (54.88, 51.85)$ and the basic reproduction number $\mathcal{R}_0 = 0.636 < 1$, which means that the epidemic spread was effectively controlled by Heilongjiang authorities.

3.3.   Contributions of different transmission routes

In this subsection, we deduce how different transmission channels from asymptomatic infection and symptomatic infection affect the epidemic spread. We also analyze the epidemic contributions of direct transmission and indirect transmission.

First, in order to separate the contribution of compartment $A$ from that of compartment $E$, we need to figure out the numbers of recovered individuals coming from asymptomatic and symptomatic individuals, respectively. Thus, we expand model (2.1) in the aim to mark the recovered individuals from different infection sources to obtain

where $H_a, \; H_e$ are the hospitalized individuals from asymptomatic compartment and exposed compartment, respectively; $R_a$ and $R_e$ are the recovered individuals from asymptomatic individuals and exposed individuals, respectively.

As above, we set parameters $(\tau, \; \delta, \; \alpha, \; \sigma, \; p)$ by (3.1) and fix the other parameters by Table 3. Numerical simulations on model (3.2), as shown in Figure 4, indicate that the infections caused by exposed individuals account for $66\%$ and the infections caused by asymptomatic individuals cover $34\%$. The result suggests that the infected individuals caused by individuals from compartment $E$ are much more than those caused by individuals from compartment $A$. Hence, it is important to adopt the contact tracing of infections from symptomatic individuals to reduce the disease spread.

Next, we infer the indirect transmission by splitting the model (2.1) to the following submodels

and

Model (3.3) consider only the transmissions from direct contacts, whereas model (3.4) includes only the indirect transmissions from environmental virus. Using the parameters set by (3.1) and Table 3, we find that the proportion of direct transmission accounts for $54\%$ and the other infections through indirect transmission covers $46\%$, which is shown in Figure 5. Thus, the indirect transmissions alone could contribute almost the half of total infections.

3.4.   Peak value estimation

It is important to study the peak value of epidemic since it provides an estimation for hospitalizing. Though $\mathcal{R}_0 < 1$, which means that the disease index $Y$ does not have a peak for $t > 0$, the numerical simulations in subsection 3.2 demonstrate that some infection compartment exhibits a peak in the early stage of epidemic, which depends on initial value. We investigate below how the peak values of epidemic disease are affected by the intensities of detection and environmental disinfection.

Let $k$ denote the coefficient of detection intensity such that the hospitalizing coefficients are described by $k\alpha$ and $k\sigma$. By setting $k = 80\%, \; 120\%, \; 150\%$ respectively and keeping the other parameter values by (3.1) and Table 3, we get the basic reproduction numbers $\mathcal{R}_0 = 0.75, \; 0.55, \; 0.45$ and peak values $A_{\rm peak} = 61.32, \; 50.44, \; 45.73$, $E_{\rm peak} = 63.51, \; 44.07, \; 36.20$ respectively. Comparing with the original peak value, we get the relative ratios of peak values for $A$ are $11.73\%, \; -8.09\%, \; -16.67\%\;$ and for $E$ are $22.50\%, \; -15.00\%, \; -30.17\%$, which imply that the intensity of detection has the more impact on exposed individuals than on asymptomatic individuals. In Figure 6, we demonstrate the progression of asymptomatic and exposed individuals, which indicates that the peak values of $A$ and $E$ get smaller and the arriving times become earlier as $k$ increases. Consequently, enhancing detection intensity is very effective on epidemic prevention and control.

To see effects of disinfection, we set $\delta = 0.45, \; 0.8, \; 1.2\; {\rm day}^{-1}$($67.5\%, \; 120\%, \; 180\%$ of original value $\delta = 0.667$) respectively and keep the other parameters by (3.1) and Table 3, we obtain the basic reproduction number $\mathcal{R}_0 = 0.637, \; 0.635, \; 0.634$ and peak values $A_{\rm peak} = 58.21, \; 53.58, \; 51.16, \; E_{\rm peak} = 55.29, \; 50.38, \; 47.42$ respectively. Comparing with the original peak value, the peak values for $A$ are increased or decreased by $6.06\%, \; -2.37\%, \; -6.78\%$ and for $E$ are increased or decreased by $6.64\%, \; -2.91\%, \; -8.53\%$ respectively. The progressions of $A$ and $E$ are shown in Figure 7. It is clear that the environmental disinfection is effective on decreasing the peak values of $A$ and $E$.

Let us now consider the case where the control measures are so loose that $\mathcal{R}_0 > 1$ or $\mathcal{R}_c > 1$. Then Theorem 2.2 means that there exists a peak for the disease index $Y(t)$. To illustrate this, we set the transition rates $\alpha = 1/14\ {\rm days}^{-1}, \sigma = 1/6\ {\rm days}^{-1}$ and disinfection rate $\delta = 1/3\ {\rm days}^{-1}$ and keep other parameters invariant, which imply that there is no detection and disinfection. Then the control reproduction number $\mathcal{R}_c = 1.78 > 1$. The curves of asymptomatic and exposed individuals and the infection index $Y(t)$ are shown in Figure 8. We see that three curves admit peaks at around time $92$ days and peak values $A_{\rm peak} = 1.32*10^6, \; E_{\rm peak} = 2.90*10^6$, which increase by $2.55*10^4,$ and $5.28*10^4$ times than original values respectively. The result implies that the weighted sum $Y(t)$ is a good estimation for predicting the emerging time of infection peak and the measures of detection and disinfection are critical in the elimination of the disease.

3.5.   Final size estimation

Based on parameters in Table 3 and (3.1), from Eq (2.4) we get the final size $z = 1.27*10^{-5}$. Then the total infection number $N_z$ can be computed by

Note that $N_z$ doesn't include $A_0, E_0$ and contains the self-healing infections $N_h$, which are omitted in statistics. By tracing infections, we get the final infection number of self-healing $N_h = 60$. As a consequence, the total infection number $N$ is

From the real statistic data in Table 2, we know that the real infection size is $N_r = 460$. Hence, the relative error is

Let us know how the final size is affected the detection rate $(k\alpha, \; k\sigma)$ and disinfection rate $\delta$ vary. As in the last subsection, we set $k = k_1 = 80\%, \; k = k_2 = 120\%, \; k = k_3 = 150\%$ and keep the rest parameters invariant. Then the corresponding final sizes become

The total infection numbers from susceptible individuals are

and the ratios are

If we set $(\delta_1, \; \delta_2, \; \delta_3) = (0.45, \; 0.8, \; 1.2)\; {\rm day}^{-1}$ ($67.5\%, \; 120\%, \; 180\%$ of original value $\delta = 0.667$) and keep other parameters invariant, then the corresponding final sizes become

The total infection numbers from susceptible individuals are

and the ratios are

Finally, we consider three special cases where $p = p_0 = 0$ (absent of asymptomatic infection), $p = p_1 = 1$ (absent of the infection from exposed individuals) and $\xi = 0$ (absent of indirect infection). From Eqs (2.12) and (2.13), we get the corresponding final sizes

The total infection numbers are

and the ratios are

From Eq (2.14) for $\xi = 0$, we obtain the final size $z_{\xi} = 4.01*10^{-6}$ and total infection number $N_{\xi} = 148 = 32.2\% N_r$. This implies that the direct transmission is the main route for epidemic transmissions.

In summary, by the numerical analysis for model (2.1), we have found how different transmission routes contribute to the infection size and how prevention measures affect the peak value and the final size of disease infection. In addition, we have demonstrated that the infection index $Y$ is a good indicator for predicting the emerging time of infection peak.

4.   Conclusions

In this paper, we have proposed the epidemic model that incorporates the disease transmission between asymptomatic infectors and susceptible individuals, and the indirect transmission from the environmental virus. These are the main transmission channels of COVID-19 resurgences in China. The disease prevalence index $Y$ is defined, which is a good indicator to predict the peak time of epidemic spread, as illustrated by Figure 8. By mathematical analysis, we have obtained the threshold condition in terms of the control reproduction number, under which the disease index admits a peak. Furthermore, the equation to compute the final size of epidemic spread is established, where the control reproduction number plays the key role. With this equation, one can calculate the contributions of different transmission routes.

We have used the data of COVID-19 resurgence in April of 2020 from Heilongjiang in China to fit the parameters of model (2.2). Figure 3 indicates that the predictive values by the mathematical model fit the real data very well. The calculations in Subsection 3.5 show that the relative error of the predictive final size of infections to the real infection size is $2.17\%$. Furthermore, our estimations from the numerical calculations to the model reveal that asymptomatic transmissions contribute $1/3$ of total infections and the symptomatic transmissions offer $2/3$ infections. We have also inferred that the indirect infection route from environmental virus is almost as dangerous as the direct transmission route of human direct contacts. By numerical calculations, we find that the peak values of the epidemic could be increased by the order $10^4$ if the prevention measures of detection and disinfection were absent. This highlights the importance of control measures during the resurgence.

In conclusion, we have proposed the mathematical model that fits the real date of COVID-19 resurgence very well, where the human-to-human transmission and indirect transmission from the environmental coronavirus play the key roles. On the basis of this model, we have derived the threshold conditions for the prevalence of the infectious disease, and obtained the mathematical equations to determine the final size of the epidemic. These provide valuable tools to design the intervention strategies to contain the epidemic according to our treatment capacity.

Acknowledgments

The research was supported by the National Natural Science Foundation of China (No.12071381). The authors are very grateful to the anonymous reviewers for their valuable suggestions that have helped to improve this paper substantially.

Conflict of interest

The authors declare there is no conflict of interest.