Stationary distribution and probability density for a stochastic SEIR-type model of coronavirus (COVID-19) with asymptomatic carriers

In this paper, we propose a stochastic SEIR-type model with asymptomatic carriers to describe the propagation mechanism of coronavirus (COVID-19) in the population. Firstly, we show that there exists a unique global positive solution of the stochastic system with any positive initial value. Then we adopt a stochastic Lyapunov function method to establish sufficient conditions for the existence and uniqueness of an ergodic stationary distribution of positive solutions to the stochastic model. Especially, under the same conditions as the existence of a stationary distribution, we obtain the specific form of the probability density around the quasi-endemic equilibrium of the stochastic system. Finally, numerical simulations are introduced to validate the theoretical findings.


Introduction
In early December of 2019, a novel coronavirus (COVID- 19) disease was reported as a major health hazard by the World Health Organization (WHO) [1]. This coronavirus is a kind of single stranded, enveloped and positive sense virus belonging to the RNA coronaviridae family [2,3]. At present, this disease has drawn much attention of researchers from all over the world and different areas. To defeat the epidemic, many researchers have made great attempts to develop various mathematical models for the transmission dynamics of this epidemic [1,[4][5][6][7][8], these models include the SEIR-type model and SEIRS-type model. In these models, the authors assumed that after a latent period, exposed individuals can show symptoms that make them easy to identify and isolate, thus they are not able to transmit the disease. However, exposed individuals can also exhibit asymptomatic after incubation, and they can continue to infect others because asymptomatic and susceptible individuals have the same set of characteristics, there is no effective way to identify them. Therefore, the classical SEIR model or SEIRS model fails to explain the discrepancy by distinguishing between the symptomatic and asymptomatic. Taking this into account, in this paper, we first develop an SEIR-type model to describe the transmission dynamics of COVID-19 which takes the following form: (1.1) Here the infectious population is split into two degenerate infectious groups and , ( ) and ( ) denote the undetected asymptomatic and symptomatic members of the infected population at time , respectively. ( ) denotes the number of individuals who are susceptible to the disease, ( ) represents the number of exposed (in the latent period) individuals and ( ) represents the number of individuals who have been infected and then removed from the possibility of being infected again at time , respectively. All parameters are strictly positive and their descriptions are given in Table 1. Note that the variable of system (1.1) does not appear explicitly in the first four equations, which implies that the individuals in the compartment do not transmit infection. By ignoring the equation for , model (1.1) reduces to the following system: (1.2)

Q. Liu and D. Jiang
According to the next generation matrix method, the basic reproduction number 0 for system (1.2) is defined by which is used to determine whether the disease occurs or not.
In addition, the corresponding dynamic behavior of system (1.2) is as follows: • If 0 ≤ 1, system (1.2) has a disease-free equilibrium 0 ( 0 , 0, 0, 0) = ( ∕ , 0, 0, 0) and it is globally asymptotically stable (GAS) in the positive invariant set , where On the other hand, it is supposed in system (1.2) that the individuals live in a constant environment. However, some parameters involved in epidemic models are always affected by the environmental noise. There has been a growing interest in considering the environmental noise in epidemiology models since it has turned out that the stochastic model can describe biological phenomena and infectious diseases in a more exact way [8][9][10][11][12][13][14][15][16]. So far, there are different ways to introduce stochastic perturbations in the model. One of the most important ways is to assume that environmental fluctuations are of white noise type which are proportional to each variable, respectively. Given the above, in the present study it is assumed that the environmental noise is separately proportional to the compartments , , and . Then corresponding to the deterministic model (1.2), the stochastic system takes the following form: Denote by 2 (R ; R + ) the family of all nonnegative functions ( ) defined on R such that they are continuously twice differentiable in . If is a vector or matrix, we use the notation ‖ ‖ to represent its norm and its transpose is defined by . For any real numbers ( = 1, … , ), we use the notation diag( 1 , … , ) to denote a -order diagonal matrix. If is an invertible matrix, we use the notation −1 to denote its inverse matrix. If is a square matrix, its determinant is denoted by | |. Moreover, if H and J are two -dimensional symmetric matrices, we define H ⪰ J ∶ H − J is at least a semi-positive definite matrix. (1.4) In view of (1.4), we can obtain that the matrix H is also positive definite if J is a positive definite matrix. The organization of this paper is summarized as follows: In Section 2, we verify that there is a unique global positive solution to the stochastic system (1.3) with any positive initial value. In Section 3, we adopt a stochastic Lyapunov function method to establish sufficient conditions for the existence and uniqueness of an ergodic stationary distribution of positive solutions to the stochastic model (1.3). In Section 4, under the same conditions as in Theorem 3.1, we obtain the accurate expression of probability density around the quasi-endemic equilibrium of the stochastic system (1.3), which reflects the strong persistence of the disease. In Section 5, numerical simulations are given to confirm our analytical findings obtained in Sections 3 and 4. Finally, a brief conclusion is given to end this paper.

Existence and uniqueness of the global positive solution
Before studying the dynamic behavior of an epidemic model, we should ensure that the solution is global and positive. The following theorem guarantees the existence and uniqueness of the global positive solution of system (1.3) with any positive initial value. Proof. Note that all the coefficients of system (1.3) are locally Lipschitz continuous, then for any initial value ( (0), (0), (0), (0)) ∈ R 4 + there is a unique local solution ( ( ), ( ), ( ), ( )) on the interval [0, ), where is an explosion time [17]. Now we validate this solution is global, i.e., to validate = ∞ a.s. To this end, let 0 ≥ 1 be sufficiently large such that (0), (0), (0) and (0) all lie within the interval [1∕ 0 , 0 ]. For each integer ≥ 0 , we define a stopping time by [17] = inf where throughout this paper we set inf ∅ = ∞ (commonly, ∅ represents the empty set). Apparently, is increasing as → ∞. Denote by s. is true, then = ∞ a.s. and ( ( ), ( ), ( ), ( )) ∈ R 4 + a.s. for all ≥ 0. Namely, to confirm the proof we need to validate ∞ = ∞ a.s. If this statement is false, then there exists a pair of constants > 0 and ∈ (0, 1) such that Accordingly, there exists an integer 1 ≥ 0 such that For any ( , , , ) ∈ R 4 + , a nonnegative 2 -function is defined by where is a positive constant to be chosen later. It is noticed that the nonnegativity of the above function is ensured due to − 1 − ln ≥ 0 for any > 0. Applying Itô's formula to differentiate , which is shown in Appendix, we have Here is a positive constant which is independent of the variables , , and . The remainder of the proof is similar to Theorem 2.1 in the literature [18] and so we omit it here. This completes the proof.

Existence of ergodic stationary distribution
In this section, we aim to establish sufficient criteria for the existence and uniqueness of an ergodic stationary distribution of positive solutions to the stochastic system (1.3). We first present some theories about the stationary distribution (see Khasminskii [19]).
Let ( ) be a regular time-homogeneous Markov process in R described by the stochastic differential equation The diffusion matrix of the process ( ) is defined as follows Proof. The proof process is divided into two steps: the first step is to prove that the uniform elliptic condition is satisfied, and the second step is to construct a nonnegative Lyapunov function which satisfies the condition ( 2 ) of Lemma 3.1.

Probability density for system (1.3)
By Theorem 3.1, we obtain that the global positive solution ( ( ), ( ), ( ), ( )) admits a unique ergodic stationary distribution (⋅) if 0 > 1. In this section, we will obtain the explicit expression of the probability density of the distribution (⋅). Although there are many works that have used the Lyapunov equation (see (4.7) below) to find the probability distribution, such as the papers [23][24][25][26], the monograph of van Kampen [27], and other works. But most of them have used numerical methods, while here we will give an analytical solution which is very interesting since there are few works on it. Firstly, two necessary transformations of system (1.3) should be introduced.
Before introducing the corresponding probability density, we still need to introduce an important definition and two necessary lemmas.

Definition 4.1 ([28]
). The characteristic polynomial of the square matrix is defined as

then is called a Hurwitz matrix if and only if
has all negative real-part eigenvalues, i.e., where the complementary definition is = 0, > . Additionally, the corresponding necessary conditions for to be a Hurwitz matrix are as follows 29,30]). For the algebraic equation 2 0 + 0 0 + 0 0 = 0, where 0 = diag(1, 0, 0, 0) and 0 is a real symmetric matrix, and the standard matrix Here 0 in this form is called the standard 1 matrix.

Probability density of stationary distribution (⋅)
Theorem 4.1. If 0 > 1, then the stationary solution ( ( ), ( ), ( ), ( )) of system (1.3) around * follows a unique log-normal probability density ( , , , ), which takes the form where the covariance matrix is positive definite, and the specific form of is given as follows.
• If 2 ≠ 0, 4 ≠ 0 and 6 ≠ 0, then ( 7 6 5 4 ) 13       Additionally, since the diffusion matrix is a constant matrix, then the probability density ( ) can be described by a Gaussian distribution through the study of Roozen [32], that is, where is a real symmetric matrix and is a positive constant satisfying the normalization condition ∫ R 4 ( ) = 1.
Since the matrices and 2 (see (4.10) below) are similar, in view of the similarity invariant of the characteristic polynomial in linear algebra, we obtain that ( ) is the similarity invariant of matrix . It is easy to get that 1 , 2 , 3 and 4 are also similarity invariants. Thus, there exists a unique standard 1 matrix of . Now we are in the position to give the specific form of and prove its positive definiteness. We divide the proof process into four steps.
Next, we will consider the following two conditions: ) .
Then we can transform Eq. (4.12) into the following form 2, we can conclude that the matrix̌0 is semi-positive definite whose explicit form is as followš Hence, the matrix 3 =̌2 3 ( 14 12   ) .
Then Eq. (4.12) can be transformed into the following equivalent form 2 it follows that the matrix 1 is semi-positive definite whose explicit form is as follows Consequently, the matrix 3 = 2 1 ( 15 11 10 ) −1 1 [( 15 11 10 ) −1 ] is also semi-positive definite and there exists a positive constant 2 such that Step 4. Consider the algebraic equation .
Next, we will consider the following two conditions:  By direct calculation and the uniqueness of the standard 1 matrix, we obtain that ) .
Then we can transform Eq. (4.13) into the following equivalent form  24 32 4 . In view of Lemma 4.2, we get that the matrix̄1 is semi-positive definite whose explicit form is as follows Accordingly, the matrix 4 = 2 2 ( 20 18 32 4 . In view of Lemma 4.2, we obtain that the matrix̃1 is semi-positive definite whose explicit form is as follows Thus, the matrix 4 = 2 3 ( 21 17   Now we are in the position to prove the matrix = 1 + 2 + 3 + 4 in Eq. (4.7) is positive definite. If 1 = 0, 3 = 0 and 5 = 0, then the covariance matrix Since the parameters 1 and ( = 1, 2, 3) are positive, so the matrix is positive definite. Similarly, we can prove that in other cases, the covariance matrix is also positive definite. Hence, according to the relationship between systems (4.1) and (4.3), we get that the stationary distribution (⋅) around * follows a unique log-normal probability density ( , , , ), which takes the form , where the specific form of can be determined by the above discussion. This completes the proof.

Numerical simulations
In this section, numerical simulations are carried out to verify the theoretical results. For the stochastic system (1.3), we adopt the Milstein higher-order method developed in [33] and the discretization form of the stochastic system is given by: where the time step > 0, 2 denote the intensity of white noises, , ( = 1, 2, 3, 4; = 1, 2, … , ) are mutually independent Gaussian random variables following the distribution (0, 1). Based on the realistic parameter values in the literature, the parameter values in all the simulations are chosen from Table 2.
Next, by numerical simulations, we mainly pay attention to validate two aspects: (i) there is a unique ergodic stationary distribution if the condition 0 > 1 holds; (ii) the existence of the probability density.  2 ) ≈ 9.3965 > 1.

Conclusion
In this paper, we formulate and analyze a stochastic SEIR-type model which is used to describe the transmission dynamics of COVID-19 in the population. Firstly, we prove that system (1.3) has a unique global positive solution with any given positive initial value. Then we use a stochastic Lyapunov function method to obtain sufficient criteria for the existence and uniqueness of an ergodic stationary distribution, which is a probability distribution with some invariant properties. In particular, under the same conditions as the existence of a stationary distribution, we get the exact expression of the probability density, which is a function that describes the probability of the output value of the random variable around the quasi-endemic equilibrium of system (1.3). Mathematically, the existence of a stationary distribution implies the weak stability in stochastic sense while the existence of the probability density of system (1.3) is more in-depth and specific than that of the stationary distribution. Biologically, the existence of a stationary distribution and probability density indicates the persistence and coexistence of all individuals.
Numerically, based on the actual parameter values in the existing literature, we obtain two important results: (i) small environmental noise makes each population fluctuate very little and hence it can retain some stochastic weak stability to some extent; (ii) we get the specific expression of the probability density around the quasi-endemic equilibrium of system (1.3). On the other hand, there are still many important topics worthy of further study. For example, in this paper, we assume that the compartment acquires permanent immunity and cannot be infected by the infectious individuals. However, in fact, some people who had been infected COVID-19 may be infected once again if they are in close contact with someone who has COVID-19 and thus they will become susceptible. With that in mind, using an SEIRS-type model to describe the transmission dynamics of COVID-19 may be more proper. But for the SEIRS-type model, because of its high dimension, when we calculate the local probability density, the discussion will become more complicated. In addition, it is also interesting to study the effects of other types of random perturbations (such as nonlinear perturbations, Poisson jumps et al.) on COVID-19 models. So far as we know, there is little literature to analyze five-dimensional epidemic models with nonlinear perturbations [38][39][40] since there are many obstacles to solving the corresponding Fokker-Planck equation due to the limitations of mathematical methods. We look forward to fully addressing it in the near future. The relevant work is now underway.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
No data was used for the research described in the article.

Acknowledgments
This work is supported by the National Natural Science Foundation of China (No. 12001090) and the Jilin Provincial Science and Technology Development Plan Project, China (No. YDZJ202201ZYTS633).

Appendix. Preliminaries of SDEs
Here we give some basic theories of stochastic differential equations (see [17] for a detailed introduction). Consider a -dimensional stochastic differential equation