Spread trend of COVID-19 epidemic outbreak in China: using exponential attractor method in a spatial heterogeneous SEIQR model

1.   Introduction

Many practical problems in real life such as epidemiology ^{[1,2,3,4,5,6,7,8,9]} and biology ^{[10,11,12,13]} can be solved by mathematical modeling. Among them, epidemiology is an important research area. The AIDS, the recent SARS, MERS, Ebola virus and outbursts of diseases such as the Zika virus are events of concern and interest to many people. During the Spring Festival of 2020, COVID-19 epidemic has spread rapidly in China, centered on Wuhan. The diffusive speed of the epidemic was unexpected, and the Chinese government quickly blocked the traffic of Wuhan to control the spread of the epidemic.

Researchers have established many types of epidemic model with diffusion. Some authors specifically to discuss the global stability of the spread of diseases ^{[14,15,16,17]}, others study the traveling wave solutions about reaction-diffusion models ^[7,8,9,18].

The spread and diffusion of epidemic disease are geographically dependent. The extent of disease transmission is different under different conditions such as altitudes, humidity and temperature. Therefore, it is very necessary to consider spatial heterogeneity in epidemiological models. In 2008, Allen et al. ^[19] proposed a diffusive susceptible–infected–susceptible (SIS) model with space-dependent, studied the existence, uniqueness, stability of the disease-free equilibrium and discussed the asymptotic behavior of the unique endemic equilibrium as $d_{S}$ approaches zero. More precisely, their results implied that if the spatial environment can be modified to include low-risk sites and the movement of susceptible individuals can be restricted, then it may be possible to eliminate the infectious disease. In some subsequent work, Peng et al. ^[20,21,22] studied the asymptotic profiles, global stability, the effects of epidemic risk and population movement for the Allen's model. In 2018, Li et al. ^[23] and Tong et al. ^[24] concerned different reaction–diffusion SIS epidemic models in a similar method, they showed the global stability of the positive constant equilibria. For spatial heterogeneous systems, the authors discussed the existence of non-constant endemic equilibrium and studied the asymptotic profile of endemic equilibria which is determined if the diffusion rate of the susceptible or infected population at $0$ or $\infty$. In 2019, Song et al. ^[25] proposed an SEIRS (Susceptible-Exposed-Infected-Recovery-Susceptible) reaction-diffusion model, where the disease transmission and recovery rates are spatially heterogeneous. They investigated the asymptotic profile of the endemic equilibrium when the diffusion rate of the susceptible individuals is small. However, we found that in the above literature, the authors have not considered the global asymptotic stability of the endemic equilibrium in the spatial heterogeneous environment.

Many diseases such as typhoid, malaria, rabies, hepatitis, AIDS, etc. cannot be completely cured at one time. These diseases will recur or relapse after treatment. Relapse, as the name implies, is a disease in the initial healing or remission phase, due to certain inducements to cause the disease to recur or relapse. Therefore, the effect of spatial factors on relapse should be taken into account. In ^[26], we discussed a class of reaction-diffusion SVIR (Susceptible-Vaccinated-Infected-Recovery) model with relapse and a varying external source in spatial heterogeneous environment as follows:

We studied the long-term dynamic behavior of this model by means of global exponential attractor theory and gradient flow method, but the diffusion coefficient $d_{S}, d_{V}, d_{I}$ and $d_{R}$ are constant number.

The density of the population can also lead to the differences in the diffusion rate of disease. For example, the recent spread of the novel coronavirus pneumonia epidemic is closely related to population density and large-scale population movement during the Spring Festival. Hence, the diffusion coefficient should be dependent on the spatial position. In this article we focus on the effects of spatial heterogeneity on diffusion. We will apply the exponential attractor method, operator semigroup theory and combine the gradient flow method to study the long-term dynamic behavior of the dissipative evolution equation. Exponential attractors can be constructed for dissipative systems which possesses a certain kind of smoothing property. And the authors in ^[27] discuss the robustness of exponentially k-dissipative dynamical systems with perturbed parameters. In this manuscript, we construct an SEIQR (Susceptible-Exposed-Infected-Quarantined-Recovery) epidemic reaction-diffusion model with relapse and a varying external source in spatially heterogeneous environment. All coefficients of this model are spatially heterogeneous, including the diffusion coefficient. Since this model is difficult to construct Lyapunov function, we will introduce the global attractor method, and give the existence condition and existence theorem of the global exponential attractor for the general dissipative evolution equation. Apply this existence theorem of the global exponential attractor and gradient flow method, the global stability of the disease-free equilibrium of the model and the persistence of the endemic equilibrium are discussed. Due to the role of spatial heterogeneity, especially the effect of relapse in spatial heterogeneity, the usual basic reproduction number $R_{0}$ has not been able to fully describe the spread of disease, but rather needs to be described by a new threshold such as eigenvalues. Finally, with the daily data on the COVID-19 published by the Chinese Health and Welfare Commission, we find that the actual spread and our theoretical results are moderate by conducting a large number of numerical simulations. Through numerical simulation, we find that the most effective way to control the spread of the epidemic is to control the contact between susceptible persons and exposed patients. Limiting contact between susceptible and infected people can only slightly reduce the number of infected people and delay the onset of time.

The organization of this paper is as follows. In Section 2, we present the existence theorem of global exponential attractor for the general dissipative evolution equation, it plays an important role in studying the global exponential attractor method of epidemic model in spatially heterogeneous environment. In Section 3, we first construct a reaction-diffusion SEIQR model with relapse and spatial heterogeneous under some assumption and prove the existence of the positive solutions. Secondly, we give the uniform boundedness and the existence of global solutions to our system. Third, we discuss the global stability and the persistence of the infection disease by using the existence theorem of global exponential attractor in Section 2. That is, we prove that the non-constant disease-free equilibrium is globally asymptotically stable when $\lambda ^{\ast } < 0$ and the epidemic disease is persisting uniformly when $\lambda ^{\ast } > 0$, where $\lambda ^{\ast }$ is a threshold. In Section 4, we simulate the spread trend and the impact on the spatial heterogeneous environment of the COVID-19 epidemic by the official data. In Section 5, we give our conclusions and some discussions.

2.   Global exponential attractor for the general dissipative evolution equation

To study the global asymptotic stability of complex models which Lyapunov functions are difficult to construct, in this Section, we will introduce the global exponential attractor method. Hence, before discussing the global dynamics of the reaction-diffusion SEIQR model with relapse and spatial heterogeneous in detail, we need to recall some concept and existing results of exponential attractor.

Definition 2.1. ([27,Definition 1.1]) Let $\mathbf{X}$ be a Banach space, $\left\{ \mathcal{Q}\left(t\right) \right\} _{t\geq 0}$ be a continuous semigroup on $\mathbf{X}$. A compact set $\mathcal{A}^{\ast }\subset \mathbf{X}$ is a global exponential attractor for $\mathcal{Q} \left(t\right)$ if

(1) it has finite fractal dimension,

(2) it is positive invariant, i.e. for arbitrary $t\geq 0$, $\mathcal{ Q}\left(t\right)\mathcal{A}^{\ast } \subset \mathcal{A}^{\ast }$,

(3) it attracts exponentially any bounded subsets of $\mathbf{X}$, i.e., for any bounded subset $\mathbf{B}\subset \mathbf{X}$, there exist constant $C\left(\left\Vert \mathbf{B}\right\Vert _{\mathbf{X}}\right) > 0$ and $\alpha > 0$, such that

If $\mathcal{A}^{\ast }$ only satisfies condition (2) and condition (3), we call it the global exponential attract set.

Definition 2.2. (^[28]) Let $\mathbf{X}\ $be a compact, connected subset of a Hilbert space $\mathbf{H}$. Let $\mathcal{Q}_{s}: \mathbf{X\rightarrow X}$ be a Lipschitz continous map. We call that $\mathcal{Q}_{s}$ has squeezing property in $\mathbf{X}$, if for some $\delta \ $in $(0, \frac{1}{4})$ there exists an orthogonal projection $P = P(\delta)$ of rank equal to $N_{0}(\delta)$ such that: for every $u$ and $v$ in $\mathbf{X}$, either

or

Theorem 2.3. (^[28]) Let $\mathcal{Q}_{s}:\mathbf{X\rightarrow X }$ be a Lipschitz continous map and $\mathcal{Q}_{s}$ have squeezing property in $\mathbf{X}$, then $\mathcal{Q}_{s}$ has an exponential fractal attractor $\mathcal{A}_{\ast }$, moreover

where $\alpha (\mathbf{X}) = \frac{\log K_{0}}{\log \frac{1}{\theta }}$, $\theta \in \left(2\delta, 1\right)$ and $K_{0}$ is a constant with respect to $\mathbf{X}$.

In the following, we prove the existence of global exponential attractor for a class of nonlinear evolution equations.

Let $\mathbf{H}_{1}, \mathbf{H}$ are Hilbert space, $\mathbf{H}_{1}\subset \mathbf{H}$ is dense and $\mathbf{H}_{1}\hookrightarrow \mathbf{H}$ is compact. $L:\mathbf{H}_{1}\rightarrow \mathbf{H}$ is a compact symmetrical sectorial operator and all eigenvalues of $L$ are

$G:\mathbf{H}_{1}\rightarrow \mathbf{H}$ is continuous.

Definition 2.4. ([29,Definition A.2.1]) A mapping $L+G: \mathbf{H}_{1}\rightarrow \mathbf{H}$ is called a gradient-type operator if there exists a $\mathbf{C}^{1}$ functional $F:\mathbf{H}_{1}\mathbf{ \rightarrow R}^{1}$ and for some constant $C_{0} > 0$ such that $DF\left(u\right) \subset \mathbf{H}$ for every $u\in \mathbf{H}_{\beta }$ for some $\beta \leq 1$, and

In this case, $F$ is called the energy functional of $L+G$, where $\mathbf{H} _{\beta }$ is the fractional power subspace generated by the sectorial operator $L$.

Lemma 2.5. Assume that condition (2.2) holds, and there exists $C_{1}, C_{2} > 0$, such that $\forall \varphi \in \mathbf{H}_{1}$, the solution $u = u(t, \varphi)$ of system (2.1) satisfies condition

where $\mathbf{H}_{\frac{1}{2}}$ is the fractional dimension subspace generated by the sectorial operator $L$. $G\left(u\right)$ is locally Lipschitz with respect to $u$. Then the system (2.1) has a global solution and exists a invariant set which exponential attracts any bounded set under the $\mathbf{H}$-norm.

Proof. By the assumption of Lemma 2.5 and $G$ is locally Lipschitz with respect to $u$, if $\forall \varphi \in \mathbf{H}_{1}$, it follows from [30,Theorem 11.3.5] that the system (2.1) has a global solution. It means that system (2.1) generates a operator semigroup $\mathcal{Q}\left(t\right)$ and$\ \mathcal{Q}\left(t\right) \varphi = u\left(t, \varphi \right)$.

We prove that for any bounded set $\mathbf{B}\subset \mathbf{H}_{1}$, there exists a invariant set $\mathbf{B}_{R}\subset \mathbf{H}$ for the operator semigroup $\mathcal{Q}\left(t\right)$. It follows from the inner product of Eq (2.1) in $\mathbf{H}$ with $u$ and by using the condition (2.3), we have that

since $\mathbf{H}_{\frac{1}{2}}\hookrightarrow \mathbf{H}$, so there exists $C_{3} > 0$, such that

Hence,

By Gronwall's inequality in differential form [26,Lemma 2.2], we obtain that

where $\alpha = 2C_{1}C_{3}^{2}$. It implies if $R^{2} > \frac{2C_{2}}{\alpha }, \varphi \in \mathbf{B}_{R}$, then

Therefore, for any $t\geq 0$, there has $u\left(t, \varphi \right) \subset \mathbf{B}_{R}$. Then, we have that $\mathbf{B}_{R}\subset \mathbf{H}$ is the invariant set.

Lemma 2.6. Assume $\mathbf{H}$ is a Hilbert space, $L:\mathbf{H} _{1}\rightarrow \mathbf{H}$ is a injective and compact symmetrical sectorial operator and $G\left(u\right)$ is Lipschitz with respect to $u$, there exists a time $t_{\ast } > 0$, such that the operator $\mathcal{Q}_{\ast } = \mathcal{Q}(t_{\ast })$ satisfies the squeezing property on $\delta < \frac{1 }{8}$ in $\mathbf{B}_{R}$, and there exists a exponential attractor $\mathcal{A}_{\ast }$ satisfies $d_{F}(\mathcal{A}_{\ast })\leq \max \{\alpha (\mathbf{B}_{R}), N_{0}\}$.

Proof. Since $L:\mathbf{H}_{1}\rightarrow \mathbf{H}$ is a injective and compact symmetrical sectorial operator, so the eigenvectors $\left\{ e_{j}\right\} _{j\in \mathbb{N}}$ corresponding to the eigenvalues $\left\{ \lambda _{j}\right\} _{j\in \mathbb{N}}$ are the complete orthonormal basis of $\mathbf{H}$, that is for any $u, v\in \mathbf{H}$, $u, v$ can be represented as follows

Then, for any ${\widetilde{\theta }}\in \left(0, \frac{1}{2}\right]$, there exists a constant $c > 0$, such that $c\left\vert \lambda _{1}\right\vert ^{ \widetilde{\theta }} > 1$, hence,

This imples that

Since $\mathbf{H}_{1}\subset \mathbf{H}$ is dense and $\mathbf{H} _{1}\hookrightarrow \mathbf{H}$ is compact, according to the properties of the fractional dimension space, $\mathbf{B}_{R}$ is compact in $\mathbf{H}$. Combine Lemma 2.5, conditions (2.4) and applying [28,Proposition 3.1] on $\mathbf{B}_{R}$, we can get that there exists a time $t_{\ast } > 0$, such that the operator $\mathcal{Q}_{\ast } = \mathcal{Q} (t_{\ast })$ satisfies the squeezing property on $\delta < \frac{1}{8}$ in $\mathbf{B}_{R}$. And there exists a exponential attractor $\mathcal{A}_{\ast }$ satisfies $d_{F}(\mathcal{A}_{\ast })\leq \max \{\alpha (\mathbf{B} _{R}), N_{0}\}$.

Next we give the existence theorem of the global exponential attractor for the general dissipative evolution equation.

Theorem 2.7. Under the hypothesis of Lemma 2.5 and Lemma 2.5 holds, then the flow $\{\mathcal{Q}(t)\}_{t\geq 0}$ generated by (2.1) admits an global exponential attractor $\mathcal{A} ^{\ast } = \mathop \cup \limits_{0 \le t \le {t_ * }} \mathcal{Q}(t)\mathcal{A} _{\ast }$, and its fractional dimension can be estimated as

Proof. For any $T > 0$, let $F\left(t, u\right) = \mathcal{Q}\left(t\right) u:\left[0, T\right] \times \mathcal{A}_{\ast }$. We verify that $F\left(t, u\right)$ is Lipschitz continuous. Since $L$ is the infinitesimal generator of $\mathcal{T}\left(t\right) = e^{tL}$ and $\mathcal{A}_{\ast }$ is compact, we have that

This shows that $\frac{\mathcal{T}\left(t\right) \varphi -\varphi }{t} -L\varphi$ is continuous on the compact set $\left[0, T\right] \times \mathcal{A}_{\ast }$, thus, it is bounded on $\left[0, T\right] \times \mathcal{A}_{\ast }$, that is, there exists a constant $\epsilon > 0$, such that

This implies that

Since $\mathcal{A}_{\ast }$ is compact and $L$ is continuous, then there exists a constant $\widetilde{C_{1}} > 0$, such that

thereby,

For any $\left(t_{1}, \varphi \right), \left(t_{2}, \varphi \right) \in \left[0, T\right] \times \mathcal{A}_{\ast }$, without loss of generality, we assume that $t_{1} > t_{2}$, then

where $\mathcal{Q}\left(t\right) \varphi = e^{tL}\varphi +\int_{0}^{t}e^{\left(t-s\right) L}G\left(\mathcal{Q}\left(s\right) \varphi \right) ds$. Then we can get that

By the boundedness of solution of system (2.1), we know that $G\left(\mathcal{Q}\left(s\right) \varphi \right)$ is bounded. Hence, there exists a constant $\widetilde{C_{2}} > 0$, such that

Since $\mathcal{T}\left(t\right) = e^{tL}$ is a continuous semigroup, then there exists a constant $M > 0$, such that

Therefore, by using the inequality (2.5), there has

where $\widetilde{C_{3}} = \max \left\{ M\widetilde{C_{1}}, M\widetilde{C_{1}} \widetilde{C_{2}}T, M\widetilde{C_{2}}\right\}$. Similarly, we can calculate that

By Gronwall's inequality, there exists a constant $\widetilde{C_{4}} > 0$, such that

Let $C = \max \left\{ \widetilde{C_{3}}, \widetilde{C_{4}}\right\}$, for any $t_{1}, t_{2}\in \left[0, T\right]$ and $\varphi _{1}, \varphi _{2}\in \mathcal{A}^{\ast }$, we have

Therefore, $F\left(t, u\right) = \mathcal{Q}\left(t\right) u:\left[0, T \right] \times \mathcal{A}_{\ast }$ is Lipschitz continuous.

By [28,Theorem 3.1], we can obtain that the flow $\{\mathcal{Q} (t)\}_{t\geq 0}$ generated by (2.1) has an exponential attractor $\mathcal{A}^{\ast } = \mathop \cup \limits_{0 \le t \le {t_ * }} \mathcal{Q}(t) \mathcal{A}_{\ast }$, and its fractional dimension can be estimated as

3.   Long-term dynamic behavior of a reaction-diffusion SEIQR epidemic model with relapse in heterogeneous environment

In this section, we construct a reaction-diffusion SEIQR (Susceptible-Exposed-Infected-Quarantined-Recovery) model with relapse in spatially heterogeneous environment. Our model is divided into five compartments, namely susceptible individuals ($S$), exposed individuals ($E$), infected individuals ($I$), quarantined individuals ($Q$) and temporary restorers ($R$). The parameters description and transfer diagram as shown below:

From Figure 1, the following system with the initial-boundary-value conditions is constructed by:

Here, $\Omega$ is a bounded domain in $\mathbb{R}^{m}\left(m\geq 1\right)$ with smooth boundary $\partial \Omega$ (when $m > 1$), $d_{S}\left(x\right), d_{E}\left(x\right), d_{I}\left(x\right), d_{R}\left(x\right) \in \mathbf{C}^{1}\left(\Omega \right)$ are positive, continuous and uniformly bounded diffusion coefficients depend on space and $\Lambda \left(x\right), \beta _{1}\left(x\right), \beta _{2}\left(x\right), \rho _{1}\left(x\right), \rho _{2}\left(x\right), \alpha \left(x\right), \gamma \left(x\right), \phi \left(x\right), \mu \left(x\right), \eta _{1}\left(x\right)$ and $\eta _{2}\left(x\right)$ are positive Hölder continuous functions on accounting for the total recruitment scale, rates of contact, relapse, incidence, quarantined, recovery, natural death and disease-related death respectively. Neumann boundary conditions $\frac{\partial S}{\partial \mathbf{n}} = \frac{ \partial E}{\partial \mathbf{n}} = \frac{\partial I}{\partial \mathbf{n}} = \frac{\partial Q}{\partial \mathbf{n}} = \frac{\partial R}{\partial \mathbf{n}} = 0$ denotes that the change rate on the boundary of the region $\Omega$ is equal to $0$. It is straightforward to verify that $\frac{SI}{S+I}\left(\frac{SE}{S+E}\right)$ is a Lipschitz continuous function of $S$ and $I\left(E\right)$ in the open first quadrant. Therefore, we can extend it to the entire first quadrant by defining it to be zero whenever $S = 0$ or $I = 0\left(E = 0\right)$. Throughout the paper, we assume that the initial value $S_{0}, E_{0}, I_{0}, Q_{0}$ and $R_{0}$ are nonnegative continuous functions on $\overline{\Omega }$, and the number of infected individuals is positive, i.e., $\int_{\Omega }I_{0}(x)dx > 0$. Specific parameters described in Table 1.

From now on, we denote that $\mathbf{H = L}^{2}\left(\Omega \right)$, $\mathbf{H}_{1}\mathbf{ = H}_{0}^{1}\left(\Omega \right) \cap \mathbf{C} ^{2, 1}\left(\Omega \right)$, and for any given continuous function $f$ on $\overline{\Omega }$,

Theorem 3.1. For each $(S_{0}(x), E_{0}(x), I_{0}(x), Q_{0}(x), R_{0} \left(x\right))\in \mathbf{C}(\overline{\Omega })$, system (3.1) exists a positive and bounded global solution $\left( S(x, t), E(x, t), I(x, t), Q(x, t), R\left(x, t\right) \right) \in \mathbf{C} ^{2, 1}(\Omega \times \left(0, \infty \right))$.

Proof. Since $L = \left(\nabla \cdot \left(d_{S}\left(x\right) \nabla \right), \nabla \cdot \left(d_{E}\left(x\right) \nabla \right), \nabla \cdot \left(d_{I}\left(x\right) \nabla \right), 0, \nabla \cdot \left(d_{R}\left(x\right) \nabla \right) \right)$ is a symmetrical sectorial operator and all eigenvalues of $L$ are

where

be quasimonotone and satisfy the locally Lipschitz conditions, so by Lemma 2.5, we can guarantees that system (3.1) exists a global solution $\left(S(x, t), E(x, t), I(x, t), Q(x, t), R\left(x, t\right) \right) \in \mathbf{C}^{2, 1}(\Omega \times \left(0, \infty \right))$. Similar to the method in [26,Lemma 2.1 and Theorem 2.2], we can easily prove that the global solution of the system (3.1) is positive. We next consider the following total population at time $t$. Define

From system (3.1), it is easy to see that

By Gronwall's inequality in differential form [26,Lemma 2.2], we can obtain that

So $U(t)\leq \max \left\{ U(0), \frac{\Lambda ^{\ast }\left\vert \Omega \right\vert }{\mu _{\ast }}\right\}$, where

This shows that $U(t) = \int_{\Omega }\left(S+E+I+Q+R\right) dx$ is bounded. By the positivity of the solution of the system (3.1), we obtain that

We denote that $K = \max \left\{ \left\Vert S\left(x, 0\right) +E\left(x, 0\right) +I\left(x, 0\right) +Q\left(x, 0\right) +R\left(x, 0\right) \right\Vert _{L^{\infty }(\Omega)}\left\vert \Omega \right\vert, \frac{ \Lambda ^{\ast }\left\vert \Omega \right\vert }{\mu _{\ast }}\right\}$, then we know

In view of [31,Theorem 1 and Corollary 1], there exists a positive constant $K^{\ast }$ depending on $K$ such that

Thus, we can obtain that $S(x, t), E(x, t), I(x, t), Q(x, t), R\left(x, t\right)$ are uniformly bounded on $\overline{\Omega }$. This implies the global solution of system (3.1) is positive and uniformly bounded.

3.1.   Threshold dynamics of the reaction-diffusion SEIQR epidemic model with relapse in heterogeneous environment

In this section, we will analyze the qualitative behavior of system (3.1). It is clearly see that system (3.1) admits a disease-free equilibrium $E^{0}\left(x\right) = (S^{0}\left(x\right), 0, 0, 0, 0)$.

In order to further study the long-term dynamic behavior of the diffusive SEIQR epidemic model with relapse in spatially heterogeneous environment, we need to prove the existence of principal eigenvalues of system (3.1). Linearizing the second, the third, the forth and the fifth equations of (3.1) at disease-free equilibrium, we get

Let $E = e^{\lambda t}\chi \left(x\right), I = e^{\lambda t}\varphi \left(x\right), Q = e^{\lambda t}\xi \left(x\right), R = e^{\lambda t}\psi \left(x\right)$, equations (3.2) can be rewritten as

Denote $\Phi \left(x\right) = \left(\chi \left(x\right), \varphi \left(x\right), \xi \left(x\right), \psi \left(x\right) \right) ^{T}, D\left(x\right) = \left[\begin{array}{llll} d_{E}\left(x\right) & 0 & 0 & 0 \\ 0 & d_{I}\left(x\right) & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & d_{R}\left(x\right) \end{array} \right]$ and

where

and $m_{ij}\left(x\right) \geq 0, i\neq j, x\in \overline{\Omega }$. Therefore, equations (3.3) can be rewritten as

By Krein–Rutman theorem, we can obtain that there exists a real eigenvalue $\lambda ^{\ast }$ of Eq (3.4) and a corresponding eigenvector $\Phi ^{\ast }\left(x\right) = \left(\chi ^{\ast }\left(x\right), \varphi ^{\ast }\left(x\right), \xi ^{\ast }\left(x\right), \psi ^{\ast }\left(x\right) \right)$ satisfying $\Phi ^{\ast }\left(x\right) > > 0$ for all $x\in \overline{\Omega }$ in the case of Neumann boundary conditions.

3.2.   Global exponential attractor and stability analysis of system (3.1)

In this section we use the global exponential attractor theory to discuss the long-term dynamic behavior of the reaction-diffusion SEIQR epidemic model with relapse in heterogeneous environment. Denote that $\mathbf{H}^{5} = \mathbf{H\times H\times H\times H\times H}$ and $\mathbf{H}_{1}^{5} = \mathbf{H }_{1}\mathbf{\times \mathbf{H}_{1}\mathbf{\times }H}_{1}\mathbf{\times H}_{1} \mathbf{\times H}_{1}$. Note that $\mathbf{H}^{5}$ and $\mathbf{H}_{1}^{5}$ are Banach spaces equipped with norm

and

In the following, we apply Lemma 2.5 and Theorem 2.7 to our system (3.1).

Theorem 3.2. System (3.1) exists a global exponential attractor $\mathcal{A}^{\ast }$, it exponential attracts any bounded set in $\mathbf{H} ^{5}$.

Proof. For the system (3.1), we first verify the condition (2.3). Since

From Theorem 3.1, we know $\left(S(x, t), E(x, t), I(x, t), Q(x, t), R\left(x, t\right) \right)$ is uniformly bounded, hence, condition (2.3) holds. In addition, since $\frac{SI}{ S+I}\left(\frac{SE}{S+E}\right)$ is a Lipschitz continuous function of $S$ and $I\left(E\right)$ in the open first quadrant, so it is not difficult for us to verify that

If we choose $\widehat{L} = \beta _{1}^{\ast }+\beta _{2}^{\ast }+\rho _{1}^{\ast }+\rho _{2}^{\ast }+\alpha ^{\ast }+\gamma ^{\ast }+\phi ^{\ast }+\mu ^{\ast }+\eta _{1}^{\ast }+\eta _{2}^{\ast }$, and denote $u = \left(S_{1}, E_{1}, I_{1}, Q_{1}, R_{1}\right), v = \left(S_{2}, E_{2}, I_{2}, Q_{2}, R_{2}\right)$, then

Hence, Lipschitz condition is well verified. As we know that

$L = \left(\nabla \cdot \left(d_{S}\left(x\right) \nabla \right), \nabla \cdot \left(d_{E}\left(x\right) \nabla \right), \nabla \cdot \left(d_{I}\left(x\right) \nabla \right), 0, \nabla \cdot \left(d_{R}\left(x\right) \nabla \right) \right)$ is a symmetrical sectorial operator and all eigenvalues of $L$ are

therefore, by Lemma 2.5, system (3.1) exists a invariant set, it exponential attracts any bounded set in $\mathbf{H}^{5}$. By Theorem 2.7, we can obtain that system (3.1) has a global exponential attractor $\mathcal{A}^{\ast }\ $and $\dim _{F}\left(\mathcal{A} ^{\ast }\right) = d_{0} < \infty$.

After getting the global exponential attractor, we discuss the stability and persists uniformly of the epidemic disease.

Theorem 3.3. The following statements are valid.

(1) If $\lambda ^{\ast } < 0$, then

in $\mathbf{H}$, and hence, the disease-free equilibrium is globally asymptotically stable.

(2) If $\lambda ^{\ast } > 0$, then there exists a function $\gamma \left(x\right) > 0$ independent of the initial data, such that any solution $(S, E, I, Q, R)$ satisfies

for $x\in \overline{\Omega }$, and hence, the disease persists uniformly

Proof. $(1)$ Suppose $\lambda ^{\ast } < 0$. We will use the comparison principle to show that $E\left(x, t\right) \rightarrow 0, I\left(x, t\right) \rightarrow 0, Q\left(x, t\right) \rightarrow 0, R\left(x, t\right) \rightarrow 0$ as $t\rightarrow \infty$ for every $x\in \Omega$. First, we observe from the system (3.1) that

Next, let us define $\left(\widetilde{E}\left(x, t\right), \widetilde{I} \left(x, t\right), \widetilde{Q}\left(x, t\right), \widetilde{R}\left(x, t\right) \right) = \left(Me^{\lambda ^{\ast }t}\chi ^{\ast }\left(x\right), Me^{\lambda ^{\ast }t}\varphi ^{\ast }\left(x\right), Me^{\lambda ^{\ast }t}\xi ^{\ast }\left(x\right), Me^{\lambda ^{\ast }t}\psi ^{\ast }\left(x\right) \right)$ where $\lambda ^{\ast } < 0, \chi ^{\ast }\left(x\right) > > 0, \varphi ^{\ast }\left(x\right) > > 0, \xi ^{\ast }\left(x\right) > > 0, \psi ^{\ast }\left(x\right) > > 0$ are the eigenvalue and eigenvectors in equations (3.3) and $M$ is chosen so large that $E\left(x, 0\right) \leq \widetilde{E}\left(x, 0\right), I\left(x, 0\right) \leq \widetilde{I} \left(x, 0\right), Q\left(x, 0\right) \leq \widetilde{Q}\left(x, 0\right), R\left(x, 0\right) \leq \widetilde{R} \left(x, 0\right)$ for every $x\in \Omega$. It can be shown that $\left(\widetilde{E}\left(x, t\right), \widetilde{I}\left(x, t\right), \widetilde{Q}\left(x, t\right), \widetilde{R}\left(x, t\right) \right)$ satisfies

By the comparison principle [26,Lemma 2.1], $E\left(x, t\right) \leq \widetilde{E}\left(x, t\right), I\left(x, t\right) \leq \widetilde{I}\left(x, t\right), Q\left(x, t\right) \leq \widetilde{Q}\left(x, t\right), R\left(x, t\right) \leq \widetilde{R}\left(x, t\right)$ for every $x\in \Omega$ and $t\geq 0$. Since $\widetilde{E}\left(x, t\right) \rightarrow 0, \widetilde{I}\left(x, t\right) \rightarrow 0, \widetilde{Q}\left(x, t\right) \rightarrow 0, \widetilde{R}\left(x, t\right) \rightarrow 0$ as $t\rightarrow \infty$ for every $x\in \Omega$, we also have that $E\left(x, t\right) \rightarrow 0, I\left(x, t\right) \rightarrow 0, Q\left(x, t\right) \rightarrow 0, R\left(x, t\right) \rightarrow 0$ as $t\rightarrow \infty$ for every $x\in \Omega$.

Next we claim $S\left(\cdot, t\right) \rightarrow S^{0}\left(x\right)$ uniformly on as $t\rightarrow \infty$. Given any small constant $\varepsilon > 0$, there exists a large time $T > 0$ such that $0\leq E(x, t), I(x, t)\leq \varepsilon$ for all $x\in \overline{\Omega }, t\geq T$. From the first equation in system (3.1), it is observed that $S$ is a super-solution to

and a sub-solution to

Denote by $w$ and $v$ the solution of system (3.5) and system (3.6), respectively. The parabolic comparison principle gives

For system (3.5), we can verify that

this means that system (3.5) satisfies condition (2.3) for $Lw+G\left(w\right) = \nabla \cdot \left(d_{S}\left(x\right) \nabla w\right) +\Lambda \left(x\right) -\left(\beta _{1}^{\ast }+\beta _{2}^{\ast }\right) \varepsilon -\mu \left(x\right) w$. It is similar to the proof of Theorem 3.2, we can obtain that system (3.5) exists a global exponential attractor $\mathcal{A}_{w}$. In addition, system (3.5) has a variational structure, the corresponding functional of the variational structure is

where

Then

so $L+G$ is a gradient type operator. From [29,Theorem A.2.2], we can prove that

where $S_{-}^{0}\left(\varepsilon, x\right)$ is the unique positive steady state of problems (3.5). Similarly, for system (3.6), we can obtain

where $S^{0}\left(x\right)$ is the unique positive steady state of problems (3.6). Furthermore, due to the arbitrariness of $\varepsilon$, it is easily checked that

Thus, our analysis implies that the $S\left(\cdot, t\right) \rightarrow S^{0}\left(x\right)$ uniformly as $t\rightarrow \infty$.

$(2)$ Since $\lambda ^{\ast } > 0$, it is observed that the solution of

is a sub-solution of the first equation in (3.1). Similar to the proof of conclusion (1), system (3.7) is also a gradient type equation. From [29,Theorem A.2.2], we can prove that

By weak maximum principle, we know that $S_{-}^{\ast }\left(x\right) > 0$ for all $x\in \Omega$. Next, let us define $\left(E_{-}\left(x, t\right), I_{-}\left(x, t\right), Q_{-}\left(x, t\right), R_{-}\left(x, t\right) \right) = \left(\varepsilon \chi ^{\ast }\left(x\right), \varepsilon \varphi ^{\ast }\left(x\right), \varepsilon \xi ^{\ast }\left(x\right), \varepsilon \psi ^{\ast }\left(x\right) \right)$ where $\chi ^{\ast }\left(x\right) > > 0, \varphi ^{\ast }\left(x\right) > > 0, \xi ^{\ast }\left(x\right) > > 0, \psi ^{\ast }\left(x\right) > > 0$ and $\varepsilon > 0$ is a sufficiently small constant. Substituting $\varepsilon \chi ^{\ast }\left(x\right), \varepsilon \varphi ^{\ast }\left(x\right), \varepsilon \xi ^{\ast }\left(x\right), \varepsilon \psi ^{\ast }\left(x\right)$ into the second, the third, the forth and the fifth equations of system (3.1), we know

and

Therefore, $\left(\varepsilon \chi ^{\ast }\left(x\right), \varepsilon \varphi ^{\ast }\left(x\right), \varepsilon \xi ^{\ast }\left(x\right), \varepsilon \psi ^{\ast }\left(x\right) \right)$ is the sub-solution of the second, the third, the forth and the fifth equations of system (3.1). We choose

$0 < \gamma \left(x\right) < \min \left\{ S_{-}^{\ast }\left(x\right), \varepsilon \chi ^{\ast }\left(x\right), \varepsilon \varphi ^{\ast }\left(x\right), \varepsilon \xi ^{\ast }\left(x\right), \varepsilon \psi ^{\ast }\left(x\right) \right\}$, we can obtain that

for $x\in \overline{\Omega }$, then it shows that the disease persists.

The synthesis of Theorem 3.2 and Theorem 3.3 can give the following results: if $\lambda ^{\ast } > 0$ and $t\rightarrow \infty$, then the positive solution of system (3.1) is globally exponential attractive and the attraction domain is $\mathcal{A}^{\ast }$. Medically speaking, the disease persists in this situation. The epidemic is globally asymptotically stable or attracted by the global exponential attractor.

4.   Numerical simulations

The COVID-19 was discovered because of Wuhan Viral Pneumonia cases in 2019. It was named by the World Health Organization on January 12, 2020. During the Spring Festival in 2020, the epidemic broke out in China with Wuhan as the center and began to diffuse. The Chinese government quickly adopted a strong response. Since January 23, 2020, Wuhan has implemented the strategy of closing the city, and various regions have adopted strategies such as restricting travel to control the spread of the epidemic. On January 30, 2020, WHO released the COVID-19 pneumonia epidemic as a public health emergency of international concern, emphasizing that travel and trade restrictions are not recommended, and reaffirming the Chinese government's control measures. While medical workers are fighting the epidemic, epidemiologists around the world have begun to use data analysis and modeling to study the spread of disease and find the fastest control methods ^{[32,33,34,35]}.

The official website of the National Health Committee of the People's Republic of China has updated the relevant data of the COVID-19 epidemic ^[36] from January 24, 2020. The official website of the WHO can check the relevant data from January 21 to the present ^[37]. The specific data is shown in the Table 2.

Judging from the existing data and various notifications of the epidemic, the COVID-19 is not yet relapse. Based on the data in Table 2, with the help of simple data analysis and calculation, we can get some important parameters in Table 3.

According to the data in Table 3 and our system (3.1), we first simulate the spread of the COVID-19 epidemic (Figure 2).

From Figure 2, we can see that the numerical simulation basically accords with the current official data, especially the exposed individuals and the infected people. In addition, we also find that the COVID-19 is currently in a period of high incidence and extremely low cure rate. In such a period, the number of quarantines in the hospital will be large, and the sickbeds in the hospital will be very tight.

The Chinese government has tried to curb the spread of the epidemic by blocking cities and traffic. The essence is to control the exposure of the disease. If we choose $\beta_2 = 10^{-8}$ in Table 3, then we can get the image in Figure 3.

From the Figure 3 we can see that merely controlling the contact between susceptible and infected people can only delay the time of getting sick, and cannot obviously reduce the number of people infected. If we want a better control, we need to control the contact rate between the susceptible and the exposed individuals, because there are more people in the exposed period and their activities are more frequent. After several numerical simulations, we found that when $\beta_1 < 0.38$, the stability of the model's endemic equilibrium began to be destroyed. If we choose $\beta_1 = 0.35$ in Table 3, then we can get a simulated image of the disease-free equilibrium of constant coefficient model (Figure 4).

We can see from the Figure 4 that if people's travel is restricted and the contact rate between susceptible and incubation period patients is controlled, the epidemic will be effectively controlled and the disease will eventually die out. At this time, the disease-free equilibrium of the COVID-19 epidemic is globally asymptotically stable. Therefore, it is very correct for the Chinese government to restrict people's access through policies such as traffic blocking and forced quarantine, and these policies may only be implemented so smoothly in China.

Next, we discuss the dynamics of the epidemic model with a spatial heterogeneous environment and relapse, so we then simulate the stability and persistence of the spatial system (3.1). From the system (3.1) we can see that all the parameters are $x-related$ functions, so we choose different functions will directly lead to different stability results. As of the end of February, a small number of relapse cases of the COVID-19 epidemic have appeared. If we choose $\beta _{1}\left(x\right) = 0.4\left\vert \sin x\right\vert, \beta _{2}\left(x\right) = \frac{0.3}{1+x}, \rho _{1}\left(x\right) = 0.4e^{-x}, \rho _{2}\left(x\right) = \frac{0.3}{2+x }, \eta _{1}\left(x\right) = 0.021e^{-x}, \eta _{2}\left(x\right) = 0.157e^{-2x}$ and select other parameters from Table 3, then we can clearly see that the disease-free equilibrium of the spatial heterogeneity epidemic is globally asymptotically stable (Figure 5).

According to our discussion in the previous Section, the diffusion coefficient is highly dependent on space. The spread of viruses and human activities are subject to conditions such as climate and temperature. From this perspective, the diffusion coefficient should have periodic characteristics. Hence, from now on, we choose $d_{S}\left(x\right) = e^{30\sin x}, d_{E}\left(x\right) = e^{15\sin x}, d_{I}\left(x\right) = e^{5\sin x}, d_{R}\left(x\right) = e^{10\sin x}$. If we then choose $\rho _{1}\left(x\right) = 0.03e^{-x}, \rho _{2}\left(x\right) = \frac{0.02}{2+x}, \eta _{1}\left(x\right) = 0.021e^{-x}$ and select other parameters from Table 3, then we can clearly see that the spatial heterogeneity epidemic is persists uniformly (Figure 6).

In Figure 6, we find that the fluctuation of the image shows an irregular and unsmooth state. This irregular and unsmooth phenomenon is caused by the spatial heterogeneity of the diffusion coefficient. However, this irregular fluctuation is still in a controllable range. And this range is also the scope of the global exponential attractor.

If we choose $\beta _{2}\left(x\right) = 0.05\sin x, \rho _{1}\left(x\right) = 0.03e^{-x}, \rho _{2}\left(x\right) = \frac{0.02}{2+x}, \eta _{1}\left(x\right) = 0.021e^{-x}$ and select other parameters from Table 3, then we can clearly see that the spatial heterogeneity epidemic is persists uniformly (Figure 7), and the disease has obvious almost periodic characteristics.

From all above simulation we can clearly see that both the disease-free equilibrium and endemic equilibrium of the constant-coefficient model are global asymptotic stability. For spatial heterogeneous models, the stability of the model is dependent on the spatial parameters and diffusion coefficients. In addition, we also see that the image whether rising or falling is very fast in the initial stage, and then fast tend to stability. It is also confirms that the solutions of system (3.1) is global exponentially attractive.

5.   Conclusions and discussion

We first prove the existence theorem of the global exponential attractor for the general dissipative evolution system. Then we formulate a reaction-diffusion SEIQR model with relapse in the spatially heterogeneous environment and prove the positivity, uniform boundedness and existence of solution of the epidemiological model. Applying the general results of global exponential attractor for dissipative evolution system in Section 2 to the model in Section 3, we can get the existence of the attractor of the epidemic model. In addition, we prove that the disease-free equilibrium is global asymptotically stable and the endemic is persisting uniformly. Due to the influence of spatial heterogeneity relapse rate, it is no longer possible to rely on the basic reproductive number to describe the spread of epidemics. Therefore, we describe the phenomenon of disease transmission through the main eigenvalue of the system. From our proof, we can find that $\lambda ^{\ast }$ has threshold characteristics. Finally, we perform numerical simulations of the spread trend of the COVID-19 epidemic in China recently. With these authoritative data, we also simulate the diffusion of the epidemic in the spatially heterogeneous environment. Simulation results show that the best way to eliminate the disease is to control the contact rate. Choosing different spatial coefficients can produce a fundamental change of the global dynamics of the system.

Acknowledgments

This work was supported by the Natural Science Foundation of Jiangsu Province, China (Grant No. BK20190578), the Fundamental Research Funds for the Central Universities, China (Grant No. JUSRP11949), the National Natural Science Foundation of China (Grant Nos. 11671180, 11671181).

Conflict of interest

The authors declare there is no conflict of interest.