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

Authors: Zhiming Li; Zhidong Teng; Changxing Ma
Journal: Complexity
 (2020)

---

## Abstract

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 numberR0 is defined and endemic equilibrium is given. Under R0>1, quasi-stationary distribution of the stochastic process is approximated by Gaussian diffusion. Residual, sensitivity, dynamical, and diffusion 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 first and most important thing is to control input, screening, treatment, and isolation.


---

## Body

## 1. Introduction

On 11 January 2020, 41 cases of pneumonia with unknown causes were reported by Wuhan Municipal Health Commission [1]. The 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]. The 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].The 2019-nCoV outbreak has received considerable attention from domestic and international scholars [7–13]. Clinical evidence of 2019-nCoV refers to [14–16]. Chan et al. [17] indicated 2019-nCoV disease can spread from person to person by a study of a family cluster. The authors of [18–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.


## 2. Methods

Suppose that the host populationN is partitioned into three compartments: susceptible, exposed (infected but symptom-free), and infectious with symptoms. Let St, Et, and It be the numbers of susceptible, exposed, and infectious individuals in the population at time t, respectively. Et+It 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. The dynamical transfer of the SEI process is demonstrated in Figure 1.Figure 1
Transfer diagram of the SEI epidemic process.In Figure1, 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. The force of infection is β1E/N and β2I/N, where β1 and β2 are defined as the effective contact per capita in the exposed and infectious periods. Thus, the incidence rate is β1SE/N+β2SI/N. Parameter σ is a transfer rate from the exposed to infectious classes.


### 2.1. SEI Markov Process

For anyN>0, we denote the process St,Et,It:t≥0 with N by a family of Markov chains SNt,ENt,INt:t≥0 that take values in state space EN≜m=s,e,i,s,e,i=0,1,2,…,N and have transition rate matrix QN≜qNm,nm,n∈EN with qNm,n representing the rate of transition from state m to n for n≠m, and qNm,m≜−qNm, where qNm=∑n≠mqNm,n<∞. For convenience, denote qNm,n=qNm,m+l with n=m+l and l=l1,l2,l3∈L≜1,0,0,−1,0,0,−1,1,0,0,−1,0,0,−1,1,0,0,−1. From Figure 1, transition rates of the SEI process are defined by six events given in Table 1.Table 1
State transitions ofSNt,ENt,INt:t≥0 and XNt:t≥0 at time t.
EventState transitionlqNm,m+lfX1,X2,X3,lInput rates,e,i⟶s+1,e,i(1, 0, 0)θNθDeath of susceptibles,e,i⟶s−1,e,i(−1, 0, 0)μsμX1Infection of susceptibles,e,i⟶s−1,e+1,i(−1, 1, 0)β1/Nse+β2/Nsiβ1X1X2+β2X1X3Death of exposeds,e,i⟶s,e−1,i(0, −1, 0)μeμX2Transfer rates,e,i⟶s,e−1,i+1(0, −1, 1)σeσX2Death of infecteds,e,i⟶s,e,i−1(0, 0, −1)μ+μ1iμ+μ1X3Note:X1=s/N,X2=e/N,andX3=i/N.Consider all states of the SEI process conditioned on nonextinction stateINt>0 at time t, denoted by SNt,ENt,INtINt>0. Let qs,e,it be condition probabilities and given by(1)qs,e,it=PSNt,ENt,INt=s,e,iINt>0,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 EN of transient states is finite and irreducible. Thus, 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 R0 and the expected population size N in steady state. When R0 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 when R0>1.
Denote(2)XNt=X1t,X2t,X3t≜SNtN,ENtN,INtN.
It can be interpreted as the population densities of the susceptible, exposed, and infectious individuals at timet. Define a function f:D⊆ℝ+3×L⟶ℝ≜−∞,+∞ such that(3)fmN,l=1NqNm,m+l,mN∈D,l∈L,FXN=∑l≠0lfXN,l,where ℝ+=0,∞. Therefore, XNt: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(4)FXN≜F1XN,F2XN,F3XN=θ−β1X1X2−β2X1X3−μX1,β1X1X2+β2X1X3−μ+σX2,σX2−μ+μ1X3.
For a density process, Ross et al. [28] can identify a deterministic analogue. Thus, the density process XNt:t≥0 should behave more deterministically as N becomes larger. The following theorem proves this point.Theorem 1.
Consider the density processXNt:t≥0 defined in (2) and limN⟶+∞XN0=x0≜X0. Then, as N⟶+∞, XNt converges uniformly in probability over finite time intervals 0,t to a unique deterministic trajectory xt satisfying x0=x0,xs∈D for 0≤s≤t and(5)dxtdt=Fxt,where xt=x1t,x2t,x3t and F⋅ is defined in (4). That is to say, for fixed t>0 and for all ϵ>0, limN⟶∞Psups≤tXNs−xs>ϵ=0.Proof.
By definition off and Table 1, it is clear that f⋅,⋅ is continuous for the first variable and therefore bounded over a compact set K⊂ℝ+3. Then, supx∈D∑llfx,l<∞,x∈D,l∈L. On the other hand, every element of F in (4) is continuous on any compact set K. Thus, F is locally Lipschitz on D; that is, there exists a constant C such that Fx−Fy<Cx−y,x,y∈D. Based on Theorem 8.1 in the study of Krutz [29], the theorem follows.
Theorem1 shows the relationship of density process XNt:t≥0 and deterministic model (5). When N is larger, we can investigate the properties of density process XNt:t≥0 by those of xt.


### 2.2. Deterministic Proportion Model

Denotext=x1t,x2t,x3t. For model (5), it is equal to the deterministic proportion model:(6)dx1tdt=θ−β1x1tx2t−β2x1tx3t−μx1t,dx2tdt=β1x1tx2t+β2x1tx3t−μ+σx2t,dx3tdt=σx2t−μ+μ1x3t.However, model (6) has no explicit solution. Thus, we discuss the local stability of equilibrium by analyzing its characteristic equations. Model (6) has an endemic equilibrium denoted by x∗=x1∗,x2∗,x3∗, satisfying(7)θ−β1x1∗x2∗−β2x1∗x3∗−μx1∗=0,β1x1∗x2∗+β2x1∗x3∗−μ+σx2∗=0,σx2∗−μ+μ1x3∗=0.And by calculation, the endemic equilibriumx∗=x1∗,x2∗,x3∗ is obtained as follows:(8)x1∗=μ+σμ+μ1β1μ+μ1+β2σ,x2∗=θβ1μ+μ1+β2σ−μμ+σμ+μ1μ+σβ1μ+μ1+β2σ,x3∗=σθβ1μ+μ1+β2σ−μμ+σμ+μ1μ+σμ+μ1β1μ+μ1+β2σ.Define the basic reproduction number as follows:(9)R0=θβ1μ+μ1+β2σμμ+σμ+μ1.Theorem 2.
IfR0>1, then the endemic equilibrium x∗=x1∗,x2∗,x3∗ is locally asymptotically stable.Proof.
The Jacobian of model (6) at point x∗ is given by(10)Jx∗=∂Fx∗=−β1x2∗−β2x3∗−μ−β1x1∗−β2x1∗β1x2∗+β2x3∗β1x1−μ+σβ2x1∗0σ−μ+μ1.
Its characteristic equation is given as follows:(11)fλ≜detλI−Jx∗=λ+β1x2∗+β2x3∗+μβ1x1∗β2x1∗−β1x2∗−β2x3∗λ−β1x1∗+μ+σ−β2x1∗0−σλ+μ+μ1=λ3+a1λ2+a2λ+a3,where(12)a1=β1x2∗+β2x3∗−β1x1∗+3μ+μ1+σ,a2=2μ+μ1+σβ1x2∗+β2x3∗−β12μ+μ1+β2σx1∗+2μμ1+σ+μ1σ+3μ2,a3=μμ+μ1+σ+μ1σβ1x2∗+β2x3∗+μ−μβ1μ+μ1+β2σx1∗.
SinceR0>1 and xi∗>0i=1,2,3, by (8), we have(13)a1=β1x2∗+β2x3∗+β2σμ+σβ1μ+μ1+β2σ+2μ+μ1>0,a2=2μ+μ1+σβ1x2∗+β2x3∗+β1μμ+μ12+β2σμμ1+2μβ1μ+μ1+β2σ>0,a3=μμ+μ1+σ+μ1σβ1x2∗+β2x3∗>0.
DenoteA=β1x2∗+β2x3∗. Then,(14)a1a2−a3>2μ+μ12μ+μ1+σA−μμ+μ1+σ+μ1σA=3μμ+μ1+μσ+μ12>0.
According toHurwitz criterion, all roots of fλ have negative real parts. Then, the endemic equilibrium x∗ is local asymptotic stability.Remark 2.
Fort⟶+∞ and larger number N, XNt converges uniformly in probability to the endemic equilibrium x∗; that is, SNt,ENt,INt converges uniformly in probability to Nx1∗,Nx2∗,Nx3∗ in the endemic phase.


### 2.3. Diffusion Approximation

Deterministic model (6) is clearly an approximation of the density process XNt:t≥0 in the endemic phase for R0>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(15)zNt=NXNt−xt.The processzNt:t≥0 is called a scaled density process. Under certain conditions, zNt converges weakly in the space of all sample paths to a Gaussian diffusion with the initial value.Theorem 3.
Consider the processzNt:t≥0 and limN⟶+∞NXN0−X0=zexist. Then, zNt converges weakly to Gaussian diffusion zt for N⟶+∞, satisfying(16)dzt=∂FXNtztdt+GXNtdwt,with initial value z, where wt:t≥0 is three-dimensional Wiener process and(17)GXN≜∑lliljfXN,l=θ+β1X1X2+β2X1X3+μX1−β1X1X2−β2X1X30−β1X1X2−β2X1X3β1X1X2+β2X1X3+μ+σX2−σX20−σX2σX2+μ+μ1X3,for XNt=X1t,X2t,X3t∈D and l∈L.Proof.
From Theorem1, trajectory of xt is contained in some compact set for N⟶+∞ because of the local asymptotic stability of the endemic equilibrium x∗. Therefore, ∂Fx is continuous and bounded to a compact set. On the other hand, by Table 1 and definition of Gxt, equation (17) yields. Further, we know that fx,l is continuous in the first variable and bounded to a compact set. Then, for l∈L, we have supx∈D∑ll2fx,l<∞. By Theorem 8.2 given by Kurtz in [29], the result follows. □
Theorem3 reveals that the scaled process zNt:t≥0 can be approximated by Gaussian diffusion zt for larger N. From Theorem 2, we get ∂FXNt⟶∂Fx∗ and GXNt⟶Gx∗ as t⟶+∞ so that zt follows a stationary O-U process given as a solution to(18)dzt=∂Fx∗ztdt+Gx∗dwt,with the initial value z. The stationary distribution of zt is multivariate normal with mean 0 and covariance matrix Σ=σij3×3 satisfying(19)∂Fx∗Σ+Σ∂Fx∗T=−Gx∗.
From the above results, we know that the O-U processzt:t≥0 has stationary mean 0 and covariance matrix Σ in the asymptotically stable case; that is, zt=NXNt−x∗∼N0,Σ. Then, the density process XNt:t≥0 can be approximated by multivariate normal distribution Nx∗,1/NΣ for larger number N and t⟶+∞. That is to say,(20)XNt=SNtN,ENtN,INtNINt>0∼Nx∗,1NΣ,for R0>1, which is diffusion approximation of the quasi-stationarity distribution qs,e,i. Further, we have(21)SNt,ENt,INtINt>0∼NNx∗,NΣ.


## 2.1. SEI Markov Process

For anyN>0, we denote the process St,Et,It:t≥0 with N by a family of Markov chains SNt,ENt,INt:t≥0 that take values in state space EN≜m=s,e,i,s,e,i=0,1,2,…,N and have transition rate matrix QN≜qNm,nm,n∈EN with qNm,n representing the rate of transition from state m to n for n≠m, and qNm,m≜−qNm, where qNm=∑n≠mqNm,n<∞. For convenience, denote qNm,n=qNm,m+l with n=m+l and l=l1,l2,l3∈L≜1,0,0,−1,0,0,−1,1,0,0,−1,0,0,−1,1,0,0,−1. From Figure 1, transition rates of the SEI process are defined by six events given in Table 1.Table 1
State transitions ofSNt,ENt,INt:t≥0 and XNt:t≥0 at time t.
EventState transitionlqNm,m+lfX1,X2,X3,lInput rates,e,i⟶s+1,e,i(1, 0, 0)θNθDeath of susceptibles,e,i⟶s−1,e,i(−1, 0, 0)μsμX1Infection of susceptibles,e,i⟶s−1,e+1,i(−1, 1, 0)β1/Nse+β2/Nsiβ1X1X2+β2X1X3Death of exposeds,e,i⟶s,e−1,i(0, −1, 0)μeμX2Transfer rates,e,i⟶s,e−1,i+1(0, −1, 1)σeσX2Death of infecteds,e,i⟶s,e,i−1(0, 0, −1)μ+μ1iμ+μ1X3Note:X1=s/N,X2=e/N,andX3=i/N.Consider all states of the SEI process conditioned on nonextinction stateINt>0 at time t, denoted by SNt,ENt,INtINt>0. Let qs,e,it be condition probabilities and given by(1)qs,e,it=PSNt,ENt,INt=s,e,iINt>0,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 EN of transient states is finite and irreducible. Thus, 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 R0 and the expected population size N in steady state. When R0 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 when R0>1.
Denote(2)XNt=X1t,X2t,X3t≜SNtN,ENtN,INtN.
It can be interpreted as the population densities of the susceptible, exposed, and infectious individuals at timet. Define a function f:D⊆ℝ+3×L⟶ℝ≜−∞,+∞ such that(3)fmN,l=1NqNm,m+l,mN∈D,l∈L,FXN=∑l≠0lfXN,l,where ℝ+=0,∞. Therefore, XNt: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(4)FXN≜F1XN,F2XN,F3XN=θ−β1X1X2−β2X1X3−μX1,β1X1X2+β2X1X3−μ+σX2,σX2−μ+μ1X3.
For a density process, Ross et al. [28] can identify a deterministic analogue. Thus, the density process XNt:t≥0 should behave more deterministically as N becomes larger. The following theorem proves this point.Theorem 1.
Consider the density processXNt:t≥0 defined in (2) and limN⟶+∞XN0=x0≜X0. Then, as N⟶+∞, XNt converges uniformly in probability over finite time intervals 0,t to a unique deterministic trajectory xt satisfying x0=x0,xs∈D for 0≤s≤t and(5)dxtdt=Fxt,where xt=x1t,x2t,x3t and F⋅ is defined in (4). That is to say, for fixed t>0 and for all ϵ>0, limN⟶∞Psups≤tXNs−xs>ϵ=0.Proof.
By definition off and Table 1, it is clear that f⋅,⋅ is continuous for the first variable and therefore bounded over a compact set K⊂ℝ+3. Then, supx∈D∑llfx,l<∞,x∈D,l∈L. On the other hand, every element of F in (4) is continuous on any compact set K. Thus, F is locally Lipschitz on D; that is, there exists a constant C such that Fx−Fy<Cx−y,x,y∈D. Based on Theorem 8.1 in the study of Krutz [29], the theorem follows.
Theorem1 shows the relationship of density process XNt:t≥0 and deterministic model (5). When N is larger, we can investigate the properties of density process XNt:t≥0 by those of xt.


## 2.2. Deterministic Proportion Model

Denotext=x1t,x2t,x3t. For model (5), it is equal to the deterministic proportion model:(6)dx1tdt=θ−β1x1tx2t−β2x1tx3t−μx1t,dx2tdt=β1x1tx2t+β2x1tx3t−μ+σx2t,dx3tdt=σx2t−μ+μ1x3t.However, model (6) has no explicit solution. Thus, we discuss the local stability of equilibrium by analyzing its characteristic equations. Model (6) has an endemic equilibrium denoted by x∗=x1∗,x2∗,x3∗, satisfying(7)θ−β1x1∗x2∗−β2x1∗x3∗−μx1∗=0,β1x1∗x2∗+β2x1∗x3∗−μ+σx2∗=0,σx2∗−μ+μ1x3∗=0.And by calculation, the endemic equilibriumx∗=x1∗,x2∗,x3∗ is obtained as follows:(8)x1∗=μ+σμ+μ1β1μ+μ1+β2σ,x2∗=θβ1μ+μ1+β2σ−μμ+σμ+μ1μ+σβ1μ+μ1+β2σ,x3∗=σθβ1μ+μ1+β2σ−μμ+σμ+μ1μ+σμ+μ1β1μ+μ1+β2σ.Define the basic reproduction number as follows:(9)R0=θβ1μ+μ1+β2σμμ+σμ+μ1.Theorem 2.
IfR0>1, then the endemic equilibrium x∗=x1∗,x2∗,x3∗ is locally asymptotically stable.Proof.
The Jacobian of model (6) at point x∗ is given by(10)Jx∗=∂Fx∗=−β1x2∗−β2x3∗−μ−β1x1∗−β2x1∗β1x2∗+β2x3∗β1x1−μ+σβ2x1∗0σ−μ+μ1.
Its characteristic equation is given as follows:(11)fλ≜detλI−Jx∗=λ+β1x2∗+β2x3∗+μβ1x1∗β2x1∗−β1x2∗−β2x3∗λ−β1x1∗+μ+σ−β2x1∗0−σλ+μ+μ1=λ3+a1λ2+a2λ+a3,where(12)a1=β1x2∗+β2x3∗−β1x1∗+3μ+μ1+σ,a2=2μ+μ1+σβ1x2∗+β2x3∗−β12μ+μ1+β2σx1∗+2μμ1+σ+μ1σ+3μ2,a3=μμ+μ1+σ+μ1σβ1x2∗+β2x3∗+μ−μβ1μ+μ1+β2σx1∗.
SinceR0>1 and xi∗>0i=1,2,3, by (8), we have(13)a1=β1x2∗+β2x3∗+β2σμ+σβ1μ+μ1+β2σ+2μ+μ1>0,a2=2μ+μ1+σβ1x2∗+β2x3∗+β1μμ+μ12+β2σμμ1+2μβ1μ+μ1+β2σ>0,a3=μμ+μ1+σ+μ1σβ1x2∗+β2x3∗>0.
DenoteA=β1x2∗+β2x3∗. Then,(14)a1a2−a3>2μ+μ12μ+μ1+σA−μμ+μ1+σ+μ1σA=3μμ+μ1+μσ+μ12>0.
According toHurwitz criterion, all roots of fλ have negative real parts. Then, the endemic equilibrium x∗ is local asymptotic stability.Remark 2.
Fort⟶+∞ and larger number N, XNt converges uniformly in probability to the endemic equilibrium x∗; that is, SNt,ENt,INt converges uniformly in probability to Nx1∗,Nx2∗,Nx3∗ in the endemic phase.


## 2.3. Diffusion Approximation

Deterministic model (6) is clearly an approximation of the density process XNt:t≥0 in the endemic phase for R0>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(15)zNt=NXNt−xt.The processzNt:t≥0 is called a scaled density process. Under certain conditions, zNt converges weakly in the space of all sample paths to a Gaussian diffusion with the initial value.Theorem 3.
Consider the processzNt:t≥0 and limN⟶+∞NXN0−X0=zexist. Then, zNt converges weakly to Gaussian diffusion zt for N⟶+∞, satisfying(16)dzt=∂FXNtztdt+GXNtdwt,with initial value z, where wt:t≥0 is three-dimensional Wiener process and(17)GXN≜∑lliljfXN,l=θ+β1X1X2+β2X1X3+μX1−β1X1X2−β2X1X30−β1X1X2−β2X1X3β1X1X2+β2X1X3+μ+σX2−σX20−σX2σX2+μ+μ1X3,for XNt=X1t,X2t,X3t∈D and l∈L.Proof.
From Theorem1, trajectory of xt is contained in some compact set for N⟶+∞ because of the local asymptotic stability of the endemic equilibrium x∗. Therefore, ∂Fx is continuous and bounded to a compact set. On the other hand, by Table 1 and definition of Gxt, equation (17) yields. Further, we know that fx,l is continuous in the first variable and bounded to a compact set. Then, for l∈L, we have supx∈D∑ll2fx,l<∞. By Theorem 8.2 given by Kurtz in [29], the result follows. □
Theorem3 reveals that the scaled process zNt:t≥0 can be approximated by Gaussian diffusion zt for larger N. From Theorem 2, we get ∂FXNt⟶∂Fx∗ and GXNt⟶Gx∗ as t⟶+∞ so that zt follows a stationary O-U process given as a solution to(18)dzt=∂Fx∗ztdt+Gx∗dwt,with the initial value z. The stationary distribution of zt is multivariate normal with mean 0 and covariance matrix Σ=σij3×3 satisfying(19)∂Fx∗Σ+Σ∂Fx∗T=−Gx∗.
From the above results, we know that the O-U processzt:t≥0 has stationary mean 0 and covariance matrix Σ in the asymptotically stable case; that is, zt=NXNt−x∗∼N0,Σ. Then, the density process XNt:t≥0 can be approximated by multivariate normal distribution Nx∗,1/NΣ for larger number N and t⟶+∞. That is to say,(20)XNt=SNtN,ENtN,INtNINt>0∼Nx∗,1NΣ,for R0>1, which is diffusion approximation of the quasi-stationarity distribution qs,e,i. Further, we have(21)SNt,ENt,INtINt>0∼NNx∗,NΣ.


## 3. Main Results

In this section, deterministic model (6) and stochastic model (21) are applied to analyze the characteristics of 2019-nCoV epidemic in Hubei province, China. The onset and death data were collected from Hubei from 15 January to 5 February 2020 [31, 32] (see Figure 2). The unknown parameters of models are estimated by least-squares method.Figure 2
Numbers of confirmed and death cases in Hubei province, from 15 January 2020 to 5 February 2020.
(a)(b)The spread of 2019-nCoV started in December 2019, in which the whole population of Hubei was 59,170,000; that is,N=59,170,000. For model (6), xiti=1,2,3 are the population proportions of susceptible, exposed, and infectious individuals of Hubei province. From Figure 2 and [32], we have I0=41 and E0=119. Obviously, S0=N−I0−E0=59,169,840. Thus,(22)x10=S0N=0.99999,x20=E0N=2.011×10−6,x30=I0N=6.929×10−7.Letθ^,β^ii=1,2,σ^,μ^, and μ^1 be the estimated values of unknown parameters θ,βii=1,2,σ,μ, and μ1, respectively. The estimated values are obtained by least-squares method (Table 2) with the initial values xi0i=1,2,3. From Table 2, we have R0=θ^β^1μ^+μ^1+β^2σ^/μ^μ^+σ^μ^+μ^1=4.2656 and x1∗=0.1139,x2∗=0.0046,andx3∗=0.0734. Since R0>1, by Theorem 2, the endemic equilibrium x∗=x1∗,x2∗,x3∗is locally asymptotically stable. Further, it reveals that 2019-nCoV infection is still spreading and will be endemic in Hubei province without effective controls. Under model (6), Figure 3 shows the fitted values of susceptible, exposed, and infectious proportions, comparing the actual values of infectious proportions.Table 2
The estimated values of model parameters by least-squares method.
ParametersDefinitionInitial valuesEstimated valuesθInput rate into the susceptible group0.0030.0034β1Contact rate in the exposed period0.4780.2800β2Contact rate in the infectious period0.3680.2936σTransfer rate from exposed to the infectious period0.5000.5554μ1Disease-caused death rate in the infectious period0.0230.0280μNatural death rate—0.007∗∗The estimation of μ is from [33].Figure 3
(a) Fitted curves ofx1t, (b) fitted curves of x2t, and (c) true and fitted curves of x3t in Hubei province, from 15 January 2020 to 5 February 2020.
(a)(b)(c)


### 3.1. Residual Analysis

Model validation is the most important step in the model building process. Residual analysis such as the coefficient of determinationR2 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). Therefore, it is reasonable to suggest that the estimators obtained here are reliable.Figure 4
95% confidence intervals of residuals for infectious proportions.Letx^3t be the fitted value of x3t and M be the number of observational data. Define(23)R2=1−∑i=1Mx3ti−x^3ti2∑i=1Mx3ti−x¯3ti2,MSE=1M∑i=1Mx3ti−x^3ti2,where x^3ti=1/M∑i=1Mx3ti. From Figure 4 and the definitions of R2 and MSE, we have R2=0.988 and MSE=1.358×10−10. Thus, model (6) can reflect the dynamics behavior of field data used in our study.


### 3.2. Sensitivity Analysis

Basic reproduction numberR0 is an important threshold value for the spread of 2019-nCoV and closely related to parameters θ, βii=1,2, σ, μ1, and μ. Given all other parameters, Figure 5 shows that R0 increases if θ, β1, or β2 is larger, but it will decrease if σ, μ1, or μ increases.Figure 5
Effects of six parameters forR0 under all other parameters given. (a)β1=0.2800,β2=0.2936,σ=0.5554,μ=0.007,andμ1=0.0280; (b)θ=0.0034,β2=0.2936,σ=0.5554,μ=0.007,andμ1=0.0280; (c)θ=0.0034,β1=0.2800,σ=0.5554,μ=0.007,andμ1=0.0280; (d)θ=0.0034,β1=0.2800,β2=0.2936,μ=0.007,andμ1=0.0280; (e)θ=0.0034,β1=0.2800,β2=0.2936,σ=0.5554,andμ=0.007; (f)θ=0.0034,β1=0.2800,β2=0.2936,σ=0.5554,andμ1=0.0280.
(a)(b)(c)(d)(e)(f)Due to uncertainty associated with the estimation of certain parameter values, it is useful to carry out sensitivity analysis to investigate how sensitiveR0 is with respect to these parameters. The sensitivity index of R0 that depends differentially on parameter τ is defined as follows:(24)ζτ=τR0dR0dτ,where τ=θ,βii=1,2,σ,μ1,andμ, and the quotient τ/R0 is introduced to normalize the index by removing the effects of units. The 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:(25)ζθ=1,ζβ1=β1μ+μ1β2σ+β1μ+μ1,ζβ2=β2σβ2σ+β1μ+μ1,ζσ=β2σμ+σ−β1σμ+μ1−β2σ2β1μ+μ1μ+σ+β2σμ+σ,ζμ1=−β2μ1σβ1μ+μ12+β2σμ+μ1,ζμ=−β1μ+μ1+β2σμ+μ12μ+σ+β2μσμ+σμ+μ1μ+σβ1μ+μ1+β2σ.For instance, takeβ1=0.2800,β2=0.2936,σ=0.5554,μ1=0.00280, and μ=0.007, then R0=1254.588θ, where the value 1254.588 is the effect of all other parameters. Thus, ζθ=1/1254.588dR0/dθ=1. According to Table 2 and the above formulas, the sensitivity indexes of R0 with respect to those parameters are listed in Table 3. The sensitivity analysis demonstrates that θ,β2,μ1, and μ are highly sensitive to changes in the value of R0. Moreover, decreasing θ and βii=1,2 means decreasing R0, different from σ,μ1, and μ.Table 3
Sensitivity indexes ofR0 about different parameters.
ParameterSensitivity indexθ1.0000β10.0567β20.9433σ−0.0442μ1−0.7546μ−1.2011


### 3.3. Dynamical and Diffusion Analyses

In this section, we mainly investigate dynamical properties ofxt and diffusion approximation of density process XNt:t≥0 in (21) when R0>1.According to the results of Section2.2, St,Et,It=Nx1t,Nx2t,Nx3t, where xiti=1,2,3 is defined in (6). Since the endemic equilibrium x∗=0.1139,0.0046,0.0734, we have Nx∗=6739463,272182,4343078. By Theorem 2, St,Et,It⟶6739463,272182,4343078 as t⟶+∞. Figure 6 provides the dynamical behavior of St,Et,It.Figure 6
Dynamical behaviors of the susceptible, exposed, and infectious individuals.
(a)(b)(c)(d)On the other hand, from (21), the process SNt,ENt,INtINt>0∼NNx∗,NΣ, where ∂Fx∗Σ+Σ∂Fx∗T=−Gx∗. After calculation, we have(26)∂Fx∗=−β^1x2∗−β^2x3∗−μ^−β^1x1∗−β^2x1∗β^1x2∗+β^2x3∗β^1x1∗−μ^+σ^β^2x1∗0σ^−μ^+μ^1=−0.0298−0.0319−0.03340.0228−0.53050.033400.5554−0.0350.Similarly, by (17), it is easy to get(27)Gx∗=0.0068−0.00260−0.00260.0052−0.00260−0.00260.0051.Then,(28)Σ=0.2418−0.0022−0.1120−0.00220.00520.0064−0.11200.00640.1753.From (21) and R0>1,(29)SNt,ENt,INt∼N6739463,272182,4343078,14307306−130174−6627042−130174307684378688−662704237868810372501.Figure7(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.Figure 7
(a) Approximation of the quasi-stationary distribution; (b) marginal density of the susceptible; (c) marginal density of the exposed; (d) marginal density of infectious individuals.
(a)(b)(c)(d)


### 3.4. Control Strategies

Note that dynamical and diffusion properties of the deterministic model and density process are based on largerN 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 cii=1,2,3 be control variables about input, screening, and isolation for treatment in order to reduce input rate θ and contact rates βii=1,2. When these control variables are introduced into model (6), we have(30)dx1tdt=θ−c1−β1−c2x1tx2t−β2−c3x1tx3t−μx1t,dx2tdt=β1−c2x1tx2t+β2−c3x1tx3t−μ+σx2t,dx3tdt=σx2t−μ+μ1x3t.Under different control settings, we calculate the endemic equilibriumx∗ and R0 according to four cases:(I) Input control variablec1>0 and other variables c2=c3=0.(II) Screening control variablec2>0 and other variables c1=c3=0.(III) Isolation and treatment variablec3>0 and other variables c1=c2=0.(IV) All control variables are nonzero.Comparing withx∗ and R0 of these cases I0∗–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 R0. Thus, 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). The detailed procedure will not be described here.Table 4
The endemic equilibriumx∗ and R0 under different control variables cii=1,2,3.
CaseControl strategyControl variablec1Control variablec2Control variablec3x∗R0I0∗—000(0.1139, 0.0046, 0.0734)4.2656I(a)0.000200(0.1139, 0.0043, 0.0678)4.0146(b)0.000400(0.1139, 0.0039, 0.0622)3.7637(c)0.001400(0.1139, 0.0021, 0.0339)2.5092(d)0.002600(0.1139, 0.0000, 0.0001)1.0037II(a)00.03000(0.1146, 0.0046, 0.0733)4.2396(b)00.18000(0.1182, 0.0046, 0.0726)4.1101(c)00.23000(0.1194, 0.0046, 0.0723)4.0669(d)00.28000(0.1207, 0.0045, 0.0721)4.0237III(a)000.0936(0.1628, 0.0040, 0.0638)2.9828(b)000.1636(0.2400, 0.0031, 0.0485)2.0234(c)000.1936(0.3013, 0.0023, 0.0364)1.6123(d)000.2336(0.4565, 0.0004, 0.0058)1.0641IV(a)0.00010.01000.0436(0.1327, 0.0042, 0.0669)3.5518(b)0.00020.02000.0936(0.1638, 0.0037, 0.0579)2.7911(c)0.00030.08000.1936(0.3147, 0.0016, 0.0253)1.4070(d)0.00040.09400.2226(0.4284, 0.0000, 0.0000)1.0003Note: the caseI0∗ corresponds to the results of Table 2.


## 3.1. Residual Analysis

Model validation is the most important step in the model building process. Residual analysis such as the coefficient of determinationR2 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). Therefore, it is reasonable to suggest that the estimators obtained here are reliable.Figure 4
95% confidence intervals of residuals for infectious proportions.Letx^3t be the fitted value of x3t and M be the number of observational data. Define(23)R2=1−∑i=1Mx3ti−x^3ti2∑i=1Mx3ti−x¯3ti2,MSE=1M∑i=1Mx3ti−x^3ti2,where x^3ti=1/M∑i=1Mx3ti. From Figure 4 and the definitions of R2 and MSE, we have R2=0.988 and MSE=1.358×10−10. Thus, model (6) can reflect the dynamics behavior of field data used in our study.


## 3.2. Sensitivity Analysis

Basic reproduction numberR0 is an important threshold value for the spread of 2019-nCoV and closely related to parameters θ, βii=1,2, σ, μ1, and μ. Given all other parameters, Figure 5 shows that R0 increases if θ, β1, or β2 is larger, but it will decrease if σ, μ1, or μ increases.Figure 5
Effects of six parameters forR0 under all other parameters given. (a)β1=0.2800,β2=0.2936,σ=0.5554,μ=0.007,andμ1=0.0280; (b)θ=0.0034,β2=0.2936,σ=0.5554,μ=0.007,andμ1=0.0280; (c)θ=0.0034,β1=0.2800,σ=0.5554,μ=0.007,andμ1=0.0280; (d)θ=0.0034,β1=0.2800,β2=0.2936,μ=0.007,andμ1=0.0280; (e)θ=0.0034,β1=0.2800,β2=0.2936,σ=0.5554,andμ=0.007; (f)θ=0.0034,β1=0.2800,β2=0.2936,σ=0.5554,andμ1=0.0280.
(a)(b)(c)(d)(e)(f)Due to uncertainty associated with the estimation of certain parameter values, it is useful to carry out sensitivity analysis to investigate how sensitiveR0 is with respect to these parameters. The sensitivity index of R0 that depends differentially on parameter τ is defined as follows:(24)ζτ=τR0dR0dτ,where τ=θ,βii=1,2,σ,μ1,andμ, and the quotient τ/R0 is introduced to normalize the index by removing the effects of units. The 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:(25)ζθ=1,ζβ1=β1μ+μ1β2σ+β1μ+μ1,ζβ2=β2σβ2σ+β1μ+μ1,ζσ=β2σμ+σ−β1σμ+μ1−β2σ2β1μ+μ1μ+σ+β2σμ+σ,ζμ1=−β2μ1σβ1μ+μ12+β2σμ+μ1,ζμ=−β1μ+μ1+β2σμ+μ12μ+σ+β2μσμ+σμ+μ1μ+σβ1μ+μ1+β2σ.For instance, takeβ1=0.2800,β2=0.2936,σ=0.5554,μ1=0.00280, and μ=0.007, then R0=1254.588θ, where the value 1254.588 is the effect of all other parameters. Thus, ζθ=1/1254.588dR0/dθ=1. According to Table 2 and the above formulas, the sensitivity indexes of R0 with respect to those parameters are listed in Table 3. The sensitivity analysis demonstrates that θ,β2,μ1, and μ are highly sensitive to changes in the value of R0. Moreover, decreasing θ and βii=1,2 means decreasing R0, different from σ,μ1, and μ.Table 3
Sensitivity indexes ofR0 about different parameters.
ParameterSensitivity indexθ1.0000β10.0567β20.9433σ−0.0442μ1−0.7546μ−1.2011


## 3.3. Dynamical and Diffusion Analyses

In this section, we mainly investigate dynamical properties ofxt and diffusion approximation of density process XNt:t≥0 in (21) when R0>1.According to the results of Section2.2, St,Et,It=Nx1t,Nx2t,Nx3t, where xiti=1,2,3 is defined in (6). Since the endemic equilibrium x∗=0.1139,0.0046,0.0734, we have Nx∗=6739463,272182,4343078. By Theorem 2, St,Et,It⟶6739463,272182,4343078 as t⟶+∞. Figure 6 provides the dynamical behavior of St,Et,It.Figure 6
Dynamical behaviors of the susceptible, exposed, and infectious individuals.
(a)(b)(c)(d)On the other hand, from (21), the process SNt,ENt,INtINt>0∼NNx∗,NΣ, where ∂Fx∗Σ+Σ∂Fx∗T=−Gx∗. After calculation, we have(26)∂Fx∗=−β^1x2∗−β^2x3∗−μ^−β^1x1∗−β^2x1∗β^1x2∗+β^2x3∗β^1x1∗−μ^+σ^β^2x1∗0σ^−μ^+μ^1=−0.0298−0.0319−0.03340.0228−0.53050.033400.5554−0.0350.Similarly, by (17), it is easy to get(27)Gx∗=0.0068−0.00260−0.00260.0052−0.00260−0.00260.0051.Then,(28)Σ=0.2418−0.0022−0.1120−0.00220.00520.0064−0.11200.00640.1753.From (21) and R0>1,(29)SNt,ENt,INt∼N6739463,272182,4343078,14307306−130174−6627042−130174307684378688−662704237868810372501.Figure7(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.Figure 7
(a) Approximation of the quasi-stationary distribution; (b) marginal density of the susceptible; (c) marginal density of the exposed; (d) marginal density of infectious individuals.
(a)(b)(c)(d)


## 3.4. Control Strategies

Note that dynamical and diffusion properties of the deterministic model and density process are based on largerN 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 cii=1,2,3 be control variables about input, screening, and isolation for treatment in order to reduce input rate θ and contact rates βii=1,2. When these control variables are introduced into model (6), we have(30)dx1tdt=θ−c1−β1−c2x1tx2t−β2−c3x1tx3t−μx1t,dx2tdt=β1−c2x1tx2t+β2−c3x1tx3t−μ+σx2t,dx3tdt=σx2t−μ+μ1x3t.Under different control settings, we calculate the endemic equilibriumx∗ and R0 according to four cases:(I) Input control variablec1>0 and other variables c2=c3=0.(II) Screening control variablec2>0 and other variables c1=c3=0.(III) Isolation and treatment variablec3>0 and other variables c1=c2=0.(IV) All control variables are nonzero.Comparing withx∗ and R0 of these cases I0∗–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 R0. Thus, 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). The detailed procedure will not be described here.Table 4
The endemic equilibriumx∗ and R0 under different control variables cii=1,2,3.
CaseControl strategyControl variablec1Control variablec2Control variablec3x∗R0I0∗—000(0.1139, 0.0046, 0.0734)4.2656I(a)0.000200(0.1139, 0.0043, 0.0678)4.0146(b)0.000400(0.1139, 0.0039, 0.0622)3.7637(c)0.001400(0.1139, 0.0021, 0.0339)2.5092(d)0.002600(0.1139, 0.0000, 0.0001)1.0037II(a)00.03000(0.1146, 0.0046, 0.0733)4.2396(b)00.18000(0.1182, 0.0046, 0.0726)4.1101(c)00.23000(0.1194, 0.0046, 0.0723)4.0669(d)00.28000(0.1207, 0.0045, 0.0721)4.0237III(a)000.0936(0.1628, 0.0040, 0.0638)2.9828(b)000.1636(0.2400, 0.0031, 0.0485)2.0234(c)000.1936(0.3013, 0.0023, 0.0364)1.6123(d)000.2336(0.4565, 0.0004, 0.0058)1.0641IV(a)0.00010.01000.0436(0.1327, 0.0042, 0.0669)3.5518(b)0.00020.02000.0936(0.1638, 0.0037, 0.0579)2.7911(c)0.00030.08000.1936(0.3147, 0.0016, 0.0253)1.4070(d)0.00040.09400.2226(0.4284, 0.0000, 0.0000)1.0003Note: the caseI0∗ corresponds to the results of Table 2.


## 4. 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 quasi-stationary 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. Theorem1 provides the relationship between the density process and the proportion model. Then, the basic reproduction number R0 and endemic equilibrium are given by the proportion model. Theorem 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 (Theorem 3). Further, it behaves like a three-dimensional O-U process, fluctuating around the endemic equilibrium of the deterministic model under R0>1.The 2019-nCoV data of Hubei province are conducted to evaluate the performance of the proposed models. Least-squares method is used to estimate unknown parameters, and residual analysis provides measures of model quality. SinceR0 is an important value to predict the spread of 2019-nCoV epidemic, sensitivity analysis reveals that R0 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 R0=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. Three 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. The strategy has not only a relative satisfactory equilibrium but also the smallest R0. Thus, 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.The work in this article bears some limitations and concerns. For instance, the time to disease extinction for 2019-nCoV is worth pursuing. Moreover, whenR0 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.



---

*Source: 9012178-2020-07-30.xml*