2019-nCoV Transmission in Hubei Province, China: Stochastic and Deterministic Analyses

Currently, a novel coronavirus (2019-nCoV) causes an outbreak of viral pneumonia in Hubei province, China. In this paper, stochastic and deterministic models are proposed to investigate the transmission mechanism of 2019-nCoV from 15 January to 5 February 2020 in Hubei province. For the deterministic model, basic reproduction number R 0 is deﬁned and endemic equilibrium is given. Under R 0 > 1, quasi-stationary distribution of the stochastic process is approximated by Gaussian diﬀusion. Residual, sensitivity, dynamical, and diﬀusion analyses of the models are conducted. Further, control variables are introduced to the deterministic model and optimal strategies are provided. Based on empirical results, we suggest that the ﬁrst and most important thing is to control input, screening, treatment, and isolation.


Introduction
On 11 January 2020, 41 cases of pneumonia with unknown causes were reported by Wuhan Municipal Health Commission [1]. e main clinical features initially include fever, cough, shortness of breath, or chest radiographs showing invasive pneumonic infiltrates in both lungs. Some patients subsequently developed pneumonia, acute respiratory distress syndrome, kidney failure, and even death [2,3]. On 12 January 2020, the World Health Organization (WHO) termed it as 2019 novel coronavirus (2019-nCoV) [4]. e National Health Commission of China announced 2019-nCoV infected pneumonia to be included in the management of statutory infectious diseases on 20 January 2020 [5]. On 11 February 2020, the WHO gave an official name to the novel coronavirus that has killed more than 1,000 people as coronavirus disease (COVID-19) [6]. e 2019-nCoV outbreak has received considerable attention from domestic and international scholars [7][8][9][10][11][12][13]. Clinical evidence of 2019-nCoV refers to [14][15][16]. Chan et al. [17] indicated 2019-nCoV disease can spread from person to person by a study of a family cluster. e authors of [18][19][20] analyzed transmission risk, the unreported number, and basic reproduction number estimations of 2019-nCoV cases in China, respectively. Mathematical modeling can be an important tool for gaining an understanding of the spread of 2019-nCoV such as stochastic and deterministic models. For stochastic versions of SI, SIS, SIR, and SIRS models, refer to [21,22]. Concerning the deterministic models, special attention has received those using ordinary differential equations and dynamical systems in their formulation [23,24].
As of 5 February 2020, more than 28,000 cases with 2019-nCoV and 600 deaths have been reported by the National Health Commission of China [5], and there were 70% of confirmed cases and 97% of deaths in Hebei province, China. In this paper, we propose a stochastic SEI process and a deterministic model to investigate 2019-nCoV transmission in Hubei province. We introduce an SEI Markovian epidemic process to analyze a density process in Section 2.1. A deterministic proportion model is proposed in Section 2.2, including an endemic equilibrium and local stability. A Gaussian diffusion approximation is considered in Section 2.3. If a scaled density process starts close to a deterministic equilibrium, this diffusion process is an Ornstein-Uhlenbeck (O-U) process. Main results are obtained in Section 3, and discussions are given in Section 4.

Methods
Suppose that the host population N is partitioned into three compartments: susceptible, exposed (infected but symptomfree), and infectious with symptoms. Let S(t), E(t), and I(t) be the numbers of susceptible, exposed, and infectious individuals in the population at time t, respectively. E(t) + I(t) is the total infected population. In the case of 2019-nCoV, there have been some clues suggesting that exposed individuals without symptom can cause many infections [25]. After one unit time, a susceptible individual can be infected through contacting with the exposed or infectious individuals and enter the E class or is still in the S class or dies. An exposed individual may have a symptom and enter the I class or still stay in E class or die. e dynamical transfer of the SEI process is demonstrated in Figure 1.
In Figure 1, the parameter θ is the input rate in the susceptible group. Parameter μ denotes the natural death rate and μ 1 is the rate of disease-caused death. e force of infection is β 1 E/N and β 2 I/N, where β 1 and β 2 are defined as the effective contact per capita in the exposed and infectious periods. us, the incidence rate is (β 1 SE/N) + (β 2 SI/N). Parameter σ is a transfer rate from the exposed to infectious classes.
Consider all states of the SEI process conditioned on nonextinction state I N (t) > 0 at time t, denoted by (S N (t), E N (t), I N (t))|I N (t) > 0. Let q s,e,i (t) be condition probabilities and given by where s, e � 0, 1, 2, . . . , N, i � 1, 2, . . . , N. Such a distribution is called quasi-stationary distribution of the SEI process. For the SEI process, the set E N of transient states is finite and irreducible.
us, the quasi-stationary distribution exists and is unique [26].

Remark 1.
According to the results of [21,27], quasi-stationary distribution has different forms depending on the basic reproductive number R 0 and the expected population size N in steady state. When R 0 is greater than 1, the distribution can be approximated by Gaussian diffusion approximation. In this work, we mainly investigate that quasi-stationary distribution is approximately normal (2) It can be interpreted as the population densities of the susceptible, exposed, and infectious individuals at time t.
Define a function f: where R + � (0, ∞). erefore, X N (t): t ≥ 0 is a density process with transition rate f given in Table 1. On the other hand, by the definition of F, we have For a density process, Ross et al. [28] can identify a deterministic analogue.
us, the density process X N (t): t ≥ 0 should behave more deterministically as N becomes larger. e following theorem proves this point.
Proof. By definition of f and Table 1, it is clear that f(·, ·) is continuous for the first variable and therefore bounded over a compact set K ⊂ R 3 + . en, sup x∈D l |l|f(x, l) < ∞, x ∈ D, l ∈ L. On the other hand, every element of F in (4) is continuous on any compact set K. us, F is locally Lipschitz on D; that is, there exists a constant C such that |F(x) − F(y)| < C|x − y|, x, y ∈ D. Based on eorem 8.1 in the study of Krutz [29], the theorem follows. eorem 1 shows the relationship of density process X N (t): t ≥ 0 and deterministic model (5). When N is larger, we can investigate the properties of density process X N (t): t ≥ 0 by those of x(t).  Complexity However, model (6) has no explicit solution. us, we discuss the local stability of equilibrium by analyzing its characteristic equations. Model (6) has an endemic equilibrium denoted by And by calculation, the endemic equilibrium is obtained as follows: Define the basic reproduction number as follows: Proof. e Jacobian of model (6) at point x * is given by Its characteristic equation is given as follows: where Note: According to Hurwitz criterion, all roots of f(λ) have negative real parts. en, the endemic equilibrium x * is local asymptotic stability. □ Remark 2. For t ⟶ + ∞ and larger number N, X N (t) converges uniformly in probability to the endemic equilibrium x * ; that is, in the endemic phase. (6) is clearly an approximation of the density process X N (t): t ≥ 0 in the endemic phase for R 0 > 1 (see Remark 2). However, it cannot reflect the fluctuations of the process around the endemic equilibrium x * . Gaussian diffusion approximation is used to analyze quasi-stationary distribution based on the work described in [29,30]. For convenience, define

Diffusion Approximation. Deterministic model
e process z N (t): t ≥ 0 is called a scaled density process. Under certain conditions, z N (t) converges weakly in the space of all sample paths to a Gaussian diffusion with the initial value.
Proof. From eorem 1, trajectory of x(t) is contained in some compact set for N ⟶ + ∞ because of the local asymptotic stability of the endemic equilibrium x * . erefore, zF(x) is continuous and bounded to a compact set. On the other hand, by Table 1 and definition of G(x(t)), equation (17) yields. Further, we know that f(x, l) is continuous in the first variable and bounded to a compact set. en, for l ∈ L, we have sup x∈D l |l| 2 f(x, l) < ∞. By eorem 8.2 given by Kurtz in [29], the result follows.
□ eorem 3 reveals that the scaled process z N (t): t ≥ 0 can be approximated by Gaussian diffusion z(t) for larger N. From eorem 2, we get zF(X N (t)) ⟶ zF(x * ) and G(X N (t)) ⟶ G(x * ) as t ⟶ +∞ so that z(t) follows a stationary O-U process given as a solution to with the initial value z. e stationary distribution of z(t) is multivariate normal with mean 0 and covariance matrix Σ � (σ ij ) 3×3 satisfying From the above results, we know that the O-U process z(t): t ≥ 0 { } has stationary mean 0 and covariance matrix Σ in the asymptotically stable case; that is, en, the density process X N (t): t ≥ 0 can be approximated by multivariate normal distribution N(x * , (1/N)Σ) for larger number N and t ⟶ +∞. at is to say, for R 0 > 1, which is diffusion approximation of the quasistationarity distribution q s,e,i . Further, we have

Residual Analysis.
Model validation is the most important step in the model building process. Residual analysis such as the coefficient of determination R 2 and mean squared error (MSE) provides measures of model quality. Figure 4 displays an error bar plot of the confidence intervals on the residuals of true and fitted values. As shown in Figure 4, residuals appear to behave randomly and the errors are relatively small for model (6). erefore, it is reasonable to suggest that the estimators obtained here are reliable. Let x 3 (t) be the fitted value of x 3 (t) and M be the number of observational data. Define where Figure 4 and the definitions of R 2 and MSE, we have R 2 � 0.988 and MSE � 1.358 × 10 − 10 . us, model (6) can reflect the dynamics behavior of field data used in our study.

Sensitivity Analysis.
Basic reproduction number R 0 is an important threshold value for the spread of 2019-nCoV and closely related to parameters θ, β i (i � 1, 2), σ, μ 1 , and μ. Given all other parameters, Figure 5 shows that R 0 increases if θ, β 1 , or β 2 is larger, but it will decrease if σ, μ 1 , or μ increases.
Due to uncertainty associated with the estimation of certain parameter values, it is useful to carry out sensitivity analysis to investigate how sensitive R 0 is with respect to these parameters. e sensitivity index of R 0 that depends differentially on parameter τ is defined as follows: where τ � θ, β i (i � 1, 2), σ, μ 1 , and μ, and the quotient τ/R 0 is introduced to normalize the index by removing the effects of units. e sensitivity index is basically the ratio of the change in output to the change in input while all other parameters remain constant. By calculation, we get the sensitivity indexes as follows: For instance, take β 1 � 0.2800, β 2 � 0.2936, σ � 0.5554, μ 1 � 0.00280, and μ � 0.007, then R 0 � 1254.588θ, where the value 1254.588 is the effect of all other parameters.

Dynamical and Diffusion Analyses.
In this section, we mainly investigate dynamical properties of x(t) and diffusion approximation of density process X N (t): t ≥ 0 in (21) when R 0 > 1.
On the other hand, from (21), the process Similarly, by (17), it is easy to get From (21) and R 0 > 1, N (6739463, 272182, 4343078), Figure 7(a) shows the approximation of the quasi-stationary distribution with a trajectory simulated from (29). In Figures 7(b)-7(d), they provide the approximation marginal trajectories of the susceptible, exposed, and infectious individuals against time, respectively.

Control Strategies. Note that dynamical and diffusion
properties of the deterministic model and density process are based on larger N and t tending to infinity. In practice, however, some good control strategies have been applied to prevent the spread of 2019-nCoV epidemic such as input control, screening, and isolation for treatment in Hubei province. Let c i (i � 1, 2, 3) be control variables about input, screening, and isolation for treatment in order to reduce input rate θ and contact rates β i (i � 1, 2). When these control variables are introduced into model (6), we have Under different control settings, we calculate the endemic equilibrium x * and R 0 according to four cases: (I) Input control variable c 1 > 0 and other variables c 2 � c 3 � 0. Comparing with x * and R 0 of these cases I * 0 -IV in Table 4, the control strategies III (d) and IV (d) are better than other strategies. On the other hand, we observe that the strategy IV (d) not only has a relatively satisfactory equilibrium x * but also produces the smallest R 0 . us, the use of strategy IV (d) is suggested for general application. Based on Table 4, dynamical properties of deterministic model (30) and diffusion properties of its density process can be obtained, similar to model (6) and (21). e detailed procedure will not be described here.

Discussions
In this paper, we propose an SEI epidemic process and a deterministic proportion model to investigate the outbreak of 2019-nCoV epidemic from 15 January 2020 to 5 February  2020 in Hubei province, China. In order to understand the property of the SEI process in endemic phase, the quasistationary distribution is of additional importance. An approximation approach is used to analyze quasi-stationary distribution based on the deterministic model.
Firstly, by scale transformation, the SEI process is equivalent to a density process with transition rate. eorem 1 provides the relationship between the density process and the proportion model. en, the basic reproduction number R 0 and endemic equilibrium are given by the proportion model. eorem 2 reveals that the density process converges uniformly in probability to the endemic equilibrium. However, it cannot reflect the fluctuations of the density process around the endemic equilibrium. To get a more accurate approximation, we make use of the Gaussian diffusion process to approximate quasi-stationary distribution ( eorem 3). Further, it behaves like a three-dimensional O-U process, fluctuating around the endemic equilibrium of the deterministic model under R 0 > 1. e 2019-nCoV data of Hubei province are conducted to evaluate the performance of the proposed models. Leastsquares method is used to estimate unknown parameters, and residual analysis provides measures of model quality. Since R 0 is an important value to predict the spread of 2019-nCoV epidemic, sensitivity analysis reveals that R 0 will decrease with smaller input rate and contact rates. Further, the sensitivity index of each parameter is given. Based on the above-estimated values, dynamical properties of the proportion model and diffusion approximation of density process are obtained for a long time, respectively. Since R 0 � 4.2656 > 1, from 15 January 2020 to 5 February 2020, 2019-nCoV infection will continue to spread and be endemic in long term if there is no effective control strategy in Hubei province (Figures 6 and 7). In order to prevent the 2019-nCoV epidemic, some good control strategies are carried out.
ree control variables are introduced into the proportion model. Table 4 reveals that the control strategy IV (d) is an optimal design, compared with other strategies. e strategy has not only a relative satisfactory equilibrium but also the smallest R 0 .
us, input control, screening, and isolation for treatment are of vital importance to prevent the spread of 2019-nCoV epidemic in Hubei province. In the practice, it has proved that the strategy is really feasible and effective to prevent the spread of 2019-nCoV epidemic. e work in this article bears some limitations and concerns. For instance, the time to disease extinction for 2019-nCoV is worth pursuing. Moreover, when R 0 is less than 1 or equal to 1, the distribution can be approximately geometric or other distributions. For the problems, we will leave these for future consideration.

Data Availability
e data used to support the results of this study are available from the National Health Commission of the People's Republic of China.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this article.