Quantifying the presymptomatic transmission of COVID-19 in the USA

1.   Introduction

The idea of invisible living creatures as agents of disease goes back at least to the writings of Aristotle (384 BC–322 BC). In the 16th century, it developed as a theory. The AIDS epidemic, the recent SARS epidemic, MERS, Ebola virus and outbursts of diseases such as the Zika virus are events of concern and interest to many people. The mechanism of transmission of infections is now known for most diseases. Some diseases transmitted by viral agents, such as influenza, measles, rubella and chicken pox, confer immunity against reinfection, others transmitted by bacteria, such as tuberculosis, meningitis, and gonorrhea, confer no immunity against reinfection. While malaria, is transmitted not directly from human to human but by vectors, usually insects that are infected by humans transmit the disease to humans. There are also diseases such as West Nile virus which are transmitted back and forth between animals and vectors ^[1]. The first known result in mathematical epidemiology is a defense of the practice of inoculation against smallpox in 1760 by Daniel Bernouilli. The first contributions to modern mathematical epidemiology are due to P.D. En'ko between 1873 and 1894 ^[2], then many scholars pay more attention to the global stability problem ^[3,4,5,6,7], bifurcation problem ^[8,9,10,11] and the optimal control problem ^[12,13] of the disease.

Recently, reaction–diffusion equations have been used by many authors in epidemiology as well as virology. Researchers have established many types of model with diffusion such as $SIR$ model ^[14,15], $SIS$ model ^[16], Lyme disease model ^[17,18], malaria model ^[19] and so on. Some authors specifically study the global stability of the spread of diseases ^{[20,21,22,23,24,25,26,27]}; others discuss the traveling wave solutions about reaction–diffusion models ^{[28,29,30,31,32,33]}. The authors of ^[34] applied energy estimates in their proofs to discussed cholera epidemic.

In the theory of epidemiology, it has been recognized that environmental heterogeneity is a significant factor in the spread of infectious diseases. Constrained by factors such as altitude, temperature, humidity, latitude, climate, and living factors, the diffusion of epidemics in different environments is vastly different. In 2008, Allen et al. ^[35] proposed a simple diffusive $SIS$ (susceptible–infected–susceptible) model with space-dependent disease transmission rate $\beta \left(x\right)$ and recovery rate $\gamma \left(x\right)$, given by

where $\Omega$ is a bounded domain in $\mathbb{R}^{m}$ and $S+I$ is the total number of individuals at $t = 0$. They studied the existence, uniqueness, stability of the disease-free equilibrium and studied 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 ^[36] discussed the asymptotic profiles of the steady states for the epidemic reaction–diffusion model (1.1). The global stability of the steady states of (1.1) was also investigated by Peng and Liu ^[37] in 2009. Four years later, Peng and Yi ^[38] studied the effects of epidemic risk and population movement on the spatiotemporal transmission of infectious disease via the $SIS$ epidemic reaction–diffusion model (1.1). In 2011, Lou and Zhao ^[19] introduced a reaction–diffusion malaria model with incubation period in the vector population. Kuniya and Wang ^[39] gave the Lyapunov functions and global stability for a spatially diffusive $SIR$ epidemic model for only one diffusion, for systems with multiple diffusions, the author did not get the method to construct the Lyapunov function. Deng and Wu ^[40] discussed the global dynamics of a $SIS$ epidemic reaction-diffusion model, which reads as follows:

The authors only discussed the attractivity and global stability of the constant coefficient model, and constructed the Lyapunov function of the constant system. Wu and Zou ^[41] analyzed the asymptotic profiles of steady states for the model (1.2) with mass action infection mechanism. In 2017, Li et al. ^[42] studied the following $SIS$ reaction-diffusion system with a linear external source:

The authors showed the uniform bounds of solutions, established the threshold dynamics in terms of the basic reproduction number and discussed the global stability of the unique constant endemic equilibrium when spatial environment is homogeneous. In particular, the asymptotic profile of endemic equilibria is determined if the diffusion rate of the susceptible or infected population at $0$ or $\infty$. Since the Lyapunov function is difficult to construct, the authors did not discuss the global stability of endemic for heterogeneous models. In 2018, Cai et al. ^[43] investigated the dynamical outcomes of a host–parasite model incorporating both horizontal and vertical transmissions in a spatial heterogeneous environment. The work of Tong and Lei ^[44] concerns an SIS epidemic reaction–diffusion model. They also failed to construct Lyapunov functions suitable for spatial heterogeneous models. Therefore, the construction of the Lyapunov function of the spatial heterogeneous model is very difficult. In discussing the long-term dynamic behavior of disease in space, we also need to try to find other methods.

On the other hand, many diseases such as hepatitis, malaria, rabies, tuberculosis and HIV/AIDS will produce relapse after being treatment. We also know that vaccination is an important means of preventing infection and relapse of these diseases. Therefore, the effect of spatial factors on relapse and vaccination should be taken into account. Inspired by the above literature, in this manuscript, we study an $SVIR$ epidemic reaction-diffusion model with relapse and a varying external source in spatially heterogeneous environment. This model has the following significant differences and difficulties compared to the above models. First, after joining the relapse, the coupling between the equations in the system is stronger. Thus, the method to prove the positivity of the solution in the prior literature can no longer be used in the coupled model. Second, because the coefficients in the model are all spatially heterogeneous and the number of equations increases, so the Lyapunov function method, which proves the global asymptotic stability of the system, is also difficult to achieve due to the increased technical difficulty of construction. For this reason, we turn to discuss the global dynamics of the system by using the attractor theory of infinite dimensional dynamical systems. It is well-known that the global attractor is an important theory to describe the long-term dynamic behavior of dissipative infinite-dimensional dynamical systems ^{[45,46,47,48,49,50,51,52,53,54]}. Finally, we conduct a large number of numerical simulations. From the simulation we can see that spatial heterogeneity has a great influence on the long-term dynamic behavior of epidemics. Reducing the relapse rate and increasing the vaccination rate can effectively control the spread of the disease.

The organization of this paper is as follows. In Section 2, we present the model under some assumption. Then by using general comparison principle of the quasi-monotone nonlinear reaction-diffusion system and the comparison principle of parabolic equation, we prove the existence of the positive solutions. We also give the uniform boundedness and the existence of global solutions to our system. In Section 3, we prove the existence theorem of global exponential attract set for our real model and apply this result to discuss the global stability and the persistence of the infection disease. That is, by applying the property of global attractor and gradient flow method, we prove that the non-constant disease-free equilibrium is globally asymptotically stable and the endemic is persisting uniformly under some new threshold conditions. To test the validity of our theoretical results, we choose some specific epidemic disease with some more practical and more definitive official data to simulate the global stability and exponential attraction of the model. The stability and global exponential attractiveness of both constant and non-constant equilibria are simulated in Section 4. In Section 5, we give our conclusions and some discussions.

2.   Positive and uniform boundedness of solutions to a reaction-diffusion $SVIR$ epidemic model with relapse in heterogeneous environment

In this section, we consider a reaction-diffusion $SVIR$ model with relapse and spatially heterogeneous environment. Our model is divided into four compartments, namely the susceptible individuals ($S$), vaccinated individuals ($V$), infected individuals ($I$) and temporary restorers ($R$). The transfer diagram and parameters description 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}, d_{V}, d_{I}, d_{R} > 0$ are diffusion coefficients and $\alpha \left(x\right), \beta _{1}\left(x\right), \beta _{2}\left(x\right), \rho \left(x\right), \phi \left(x\right), \mu \left(x\right), \eta _{1}\left(x\right), \eta _{2}\left(x\right)$ and $\Lambda \left(x\right)$ are positive Hölder continuous functions on accounting for the rates of vaccination, disease transmission, relapse, recovery, natural death, disease-related death and total recruitment scale respectively. Neumann boundary conditions $\frac{\partial S}{\partial {\bf{n}}} = \frac{ \partial V}{\partial {\bf{n}}} = \frac{\partial I}{\partial {\bf{n}}} = \frac{\partial R}{\partial {\bf{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{VI}{V+I}\right)$ is a Lipschitz continuous function of $S$ $\left(V\right)$ and $I$ in the open first quadrant. Therefore, we can extend it to the entire first quadrant by defining it to be zero whenever $S = 0$ $\left(V = 0\right)$ or $I = 0$. Throughout the paper, we assume that the initial value $S_{0}(x), V_{0}(x), I_{0}(x)$ and $R_{0}(x)$ 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.

With the development of science and technology and the improvement of medical standards, vaccines play an important role in the prevention and control of disease transmission. So we considered a group of people who were vaccinated with preventive vaccines in the model (2.1). In the model, $I$ denotes the infected person, and $R$ denotes the person who is temporarily recovered by treatment after the infection but may also relapse. These two groups of people carry the disease-causing virus, so we need to consider their disease-induced rates. $S$ denotes a susceptible population, and $V$ denotes an uninfected person who has been vaccinated. They are not ill before they are exposed to pathogenic bacteria or infected people, so we only need to consider the natural mortality of these two types of people. There are many diseases in real life that have such properties, such as hepatitis B, smallpox, etc.,

In the following, we discuss the positivity and boundedness of the solution of system (2.1). In many articles, the positivity, boundedness, and existence of solutions of the model have been discussed, such as models that lose immunity in ^[55,56,57]. The authors in ^[55] construct a model of ordinary differential equations, which is naturally less difficult for the proof of positive, boundedness, and existence of solutions. For example, the positive problem of the solution can be easily obtained through contradiction analysis even without using the comparison principle. The model in our manuscript is diffused, and the difficulty of analysis is naturally higher than that of the ODE model. Models in ^[56] and ^[57] are diffusion models, and the model of ^[57] assumes that all diffusion coefficients are equal at the time of construction. In ^[56], when discussing the existence of the solution of the model (Theorem 2.1 in ^[56]), the condition that the diffusion coefficients are equal is added. We all know that it is very easy to add the equations in the system when the diffusion coefficients are equal, and then researchers can discuss positivity, boundedness and existence of the solution for this equation about the rate of change of the total population. However, the model diffusion coefficients in our manuscripts are different, we refuse to add any conditions to the diffusion coefficient from beginning to end. Therefore, from the perspective of the diffusion system, our model is more general. These are places where our models and discussion methods differ from other existing results. Firstly, we give an comparison theorem for semilinear parabolic differential equation that we need to use later.

Denote

where ${\bf{n}}$ is the unit outside normal vector of $\partial \Omega$ and $b_{i}\left(x\right) \geq 0$. The following lemma is similar to the Theorem 4.2.9 in Chinese books ^[58], but for the convenience of the reader, we give its conclusions and a brief proof.

Lemma 2.1. Let $u, v\in {\bf{C}}^{2, 1}\left(Q_{T}\right) \cap {\bf{C}}\left(\overline{Q}_{T}\right)$, satisfy

and

Then

(1) $u\left(x, t\right) \geq v\left(x, t\right)$ in $Q_{T}$.

(2) If there exist $\left(x_{0}, t_{0}\right) \in Q_{T}$ and $\mu :0\leq \mu \leq m$, such that $\left(u-v\right) _{\mu }\left(x_{0}, t_{0}\right) = 0$, then for any $\left(x, t\right) \in \overline{Q} _{t_{0}}$, there has $\left(u-v\right) _{\mu }\left(x, t\right) = 0$.

Proof. Let $w = u-v$, then by using differential mean value theorem for multivariate functions

In view of $\frac{\partial f_{i}}{\partial u_{j}}\geq 0, i\neq j, i, j = 1, 2, ..., m$, therefore, $h_{ij} = -\frac{\partial f_{i}}{\partial u_{j}} \leq 0$. If $x\in \overline{\Omega }$, we have $w\left(x, 0\right) \geq 0$ and $B_{i}w_{i}|_{S_{T}}\geq 0, i = 1, 2, ..., m$. Then by the strong maximum principle in the form of equation derivatives, we can prove that $w\left(x, t\right) \geq 0$, that is $u\left(x, t\right) \geq v\left(x, t\right)$ in $Q_{T}$. Finally, if there exist $\left(x_{0}, t_{0}\right) \in Q_{T}$ and $\mu :0\leq \mu \leq m$, such that $w_{\mu }\left(x_{0}, t_{0}\right) = 0$, then for any $\left(x, t\right) \in \overline{Q}_{t_{0}}$, there has $w_{\mu }\left(x, t\right) = 0$, that is $\left(u-v\right) _{\mu }\left(x, t\right) = 0$.

The above lemma can also be derived indirectly from [59,Theorem 7.3.4]. Let ${\bf{X}}: = {\bf{L}}^{2}\left(\Omega \right)$ be a Banach space, $L_{1}, L_{2}, L_{3}$ and $L_{4}$ be a linear operator on $D\left(L\right) = \left(D\left(L_{1}\right), D\left(L_{2}\right), D\left(L_{3}\right), D\left(L_{4}\right) \right) \subset {\bf{X}}$ defined by

Then, we can see that $L_{1}, L_{2}, L_{3}, L_{4}$ are respectively the infinitesimal generator of strongly continuous semigroup $\left\{ e^{tL_{1}}\right\} _{t\geq 0}, \left\{ e^{tL_{2}}\right\} _{t\geq 0}, \left\{ e^{tL_{3}}\right\} _{t\geq 0}, \left\{ e^{tL_{4}}\right\} _{t\geq 0}$ in ${\bf{H}}^{2}\left(\Omega \right)$. Let $\left({\bf{H}}^{2}\left(\Omega \right) \right) ^{4} = {\bf{H}}^{2}\left(\Omega \right) {\bf{\times H}}^{2}\left(\Omega \right) \times {\bf{H}}^{2}\left(\Omega \right) \times {\bf{H}} ^{2}\left(\Omega \right)$ and define $L:\left({\bf{H}}^{2}\left(\Omega \right) \right) ^{4}{\bf{\rightarrow }}\left({\bf{H}}^{2}\left(\Omega \right) \right) ^{4}, \forall u: = \left(S, V, I, R\right) \in \left({\bf{H}}^{2}\left(\Omega \right) \right) ^{4}$, then

is also the infinitesimal generator of a strongly continuous semigroup $\left\{ e^{tL}\right\} _{t\geq 0}$ in ${\bf{Y: = }}\left({\bf{H}}^{2}\left(\Omega \right) \right) ^{4}$, where

Note that ${\bf{Y}}$ is a Banach space equipped with norm

Let $G$ be a nonlinear operator on ${\bf{Y}}$ defined by

where

Let $F = \left(f_{1}, f_{2}, f_{3}, f_{4}\right)$, where

Using the operators $L$ and $F$ defined by above, we can rewrite the problem (2.1) into the following abstract form in ${\bf{Y}}$.

Because $L$ is the generator of a $C_{0}$-semigroup of contractions on ${\bf{Y}}$ and $F$ is locally Lipschitz with respect to $u$, apply [60,Corollary 11.1.2] and [60,Corollary 11.3.1], we know that there exists at least one saturated solution $u:\left[0, T\right) \rightarrow {\bf{Y}}$, which is a classical solution of system (2.2), such that

Moreover, if $S(x, 0), V(x, 0), I(x, 0), R(x, 0)\in {\bf{C}}^{2}(\Omega)$, then $u\left(x, t\right) \in {\bf{C}}^{2, 1}(\Omega \times \left(0, T\right])$.

In the following, we discuss the positivity of the solution of system (2.2) and the existence of the solution of system (2.1).

Theorem 2.2. Assume that $(S(x, t), V(x, t), I(x, t), R(x, t))$ be a solution of system (2.2), and $S(x, t), V(x, t), I(x, t), R(x, t)\in {\bf{C }}(\overline{\Omega }\times \left(0, T\right])\cap {\bf{C}}^{2, 1}(\Omega \times \left(0, T\right])$ where $T$ is the maximal existing time. If $S(x, 0) > 0, V(x, 0) > 0, I(x, 0) > 0, R(x, 0) > 0$ for $x\in \Omega$, then the solution $S(x, t) > 0, V(x, t), > 0, I(x, t) > 0, R(x, t) > 0, (x, t)\in \Omega \times \left(0, T\right]$ and $u\left(t\right) = (S(x, t), V(x, t), I(x, t), R(x, t))$ is a solution ofsystem (2.1).

Proof. We denote that

It is obviously that $\widetilde{f_{1}}\left(S, V, I, R\right) \leq f_{1}\left(S, V, I, R\right)$ and $0$ is a solution of equation

If $S$ is the solution of the equation

then by using [61,Theorem 2.2.1], we get $S(x, t) > 0$. Similarly, it is obviously that $\widetilde{f_{2}}\left(S, V, I, R\right) \leq f_{2}\left(S, V, I, R\right)$ and $0$ is a solution of equation

If $V$ is the solution of the equation

then by using [61,Theorem 2.2.1], we get $V(x, t) > 0$. At this time, $S$ and $V$ are just the solutions of the first and the second equations system (2.2). Then for the third and the forth equations of system (2.2), we know that

therefore, $\frac{\partial f_{3}}{\partial R}\geq 0, \frac{\partial f_{4}}{ \partial I}\geq 0$. By the comparison principle Lemma 2.1, we can obtain $I(x, t) > 0$ and $R(x, t) > 0$. Because under the condition of $S, V, I, R > 0$, we have $F = G$, therefore, $u\left(t\right)$ is also the solution of (2.1).

From now on, for any given continuous function $f$ on $\overline{\Omega }$, we denote

In the following, we need to give the Gronwall type differential inequality. Moreover, we do not make assumptions about the positivity of all the functions that appear in the lemma. This lemma will play an important role in the proof of the following lemmas and theorems.

Lemma 2.3. $u\left(t\right), c\left(t\right), h\left(t\right)$ are three real functions defined on $\left[a, T\right)$, $u\left(t\right)$ is absolutely continuous on $\left[a, T\right)$, $c\left(t\right)$ is integrable on $\left[a, T\right)$. If the following inequality holds,

then

Theorem 2.4. Assume that $(S(x, t), V(x, t), I(x, t), R(x, t))$ be a solution of system (2.1), and $S(x, t), V(x, t), I(x, t), R(x, t)\in {\bf{C }}(\overline{\Omega }\times \left(0, T\right])\cap {\bf{C}}^{2, 1}(\Omega \times \left(0, T\right])$ where $T$ is the maximal existing time. If $S(x, 0) > 0, V(x, 0) > 0, I(x, 0) > 0, R(x, 0) > 0$ for $x\in \Omega$, then the solution $\left(S(x, t), V(x, t), I(x, t), R(x, t)\right)$ uniformly bounded on $\overline{ \Omega }$.

Proof. We consider the following total population at time $t$. Define

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

By Lemma 2.3, 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+V+I+R\right) dx$ is bounded. By Theorem 2.2, we obtain that

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

In view of [62,Lemma 2.1], there exists a positive constant $K^{\ast }$ depending on $K$ such that

Thus, we can obtain that $S(x, t), V(x, t), I(x, t), R(x, t)$ are uniformly bounded on $\overline{\Omega }$. This implies the boundedness of the solution of system (2.1).

By [60,Corollary 11.3.2] we can deduce that the saturated solution $u\left(t\right) = \left(S(x, t), V(x, t), I(x, t), R(x, t)\right)$ is either global or $\underset{ t\rightarrow +\infty }{\lim }\left\Vert u(t)\right\Vert _{\bf{Y} } = +\infty$. From Theorem 2.4 we know that $u\left(t\right)$ is bounded, hence, $u\left(t\right)$ is a global solution.

3.   Threshold dynamics of the reaction-diffusion $SVIR$ epidemic model with relapse in heterogeneous environment

In this section, we will analyze the qualitative behavior of system (2.1). It is clearly seen that system (2.1) admits a disease-free equilibrium $E^{0}\left(x\right) = (S^{0}\left(x\right), V^{0}\left(x\right), 0, 0)$, then we prove that the existence of principal eigenvalues of system (2.1). Linearizing the second and the third equations of (2.1) at disease-free equilibrium, we get

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

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

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

By [59,Theorem 7.6.1], we can obtain that there exists a real eigenvalue $\lambda ^{\ast }$ of Eq (3.3) and a corresponding eigenvector $\Phi ^{\ast }\left(x\right)$ satisfying $\Phi ^{\ast }\left(x\right) > > 0$ for all $x\in \overline{\Omega }$ in the case of Neumann boundary conditions.

3.1.   Long-term dynamics of the reaction-diffusion $SVIR$ epidemic model

Before discussing the global dynamics of the reaction-diffusion $SVIR$ model with relapse and spatial heterogeneous in detail, we need to recall some concept and existing results of exponential attractor and measure of noncompactness.

Let ${\bf{X}}$ be a Banach space with the decomposition

and denote orthogonal projector by $P:{\bf{X}}\rightarrow {\bf{X}}_{1}$ and $(I-P):{\bf{X}}\rightarrow {\bf{X}}_{2}$. In addition, let $\left\{ Q\left(t\right) \right\} _{t\geq 0}$ be a continuous semigroup on ${\bf{X }}$. The following Condition (C*) introduced in ^[53] is an important condition for verifying the global exponential attraction.

Condition (C*): For any bounded set ${\bf{B}}\subset {\bf{X}}$ there exist positive constants $t_{\bf{B}}, C$ and $\alpha$ such that for any $\varepsilon > 0$ there exists a finite dimensional subspace ${\bf{ X}}_{1}\subset {\bf{X}}$ satisfies

where $P:{\bf{X}}\rightarrow {\bf{X}}_{1}$ is a orthogonal projector.

In the following, we denote that ${\bf{H = L}}^{2}\left(\Omega \right) \cap {\bf{C}}^{2, 1}\left(\Omega \right)$, ${\bf{H}}_{1}{\bf{ = H}} _{0}^{1}\left(\Omega \right) \cap {\bf{C}}^{2, 1}\left(\Omega \right), {\bf{H}}^{4} = {\bf{H\times H\times H\times H}}$ and ${\bf{H}}_{1}^{4} = {\bf{H}}_{1}{\bf{\times H}}_{1}{\bf{\times H}}_{1}{\bf{\times H}} _{1}$. Note that ${\bf{H}}^{4}$ and ${\bf{H}}_{1}^{4}$ are Banach spaces equipped with norm

and

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

In this case, $F$ is called the energy functional of $L+G$.

Theorem 3.2. System (2.1) exists a global exponential attract set $\mathcal{A}^{\ast }$, it exponential attracts any bounded set in ${\bf{H}}^{4}$.

Proof. It is known by the existence of the global solution that if $\forall \varphi = u\left(0\right) = \left(S_{0}\left(x\right), V_{0}\left(x\right), I_{0}\left(x\right), R_{0}\left(x\right) \right) ^{T}\in {\bf{H}}_{1}^{4}$, then system (2.1) has a global solution $u = \left(S, V, I, R\right) ^{T} \in {\bf{C}}^{0}\left(\left[0, \infty \right), {\bf{H}}^{4}\right)$. It means that system (2.1) generates a operator semigroup $Q\left(t\right) = \left(Q_{1}\left(t\right), Q_{2}\left(t\right), Q_{3}\left(t\right), Q_{4}\left(t\right) \right) ^{T}$ and$\ Q\left(t\right) \varphi = u\left(t, \varphi \right)$.

We first prove that the operator semigroup $Q\left(t\right)$ has an absorbing set ${\bf{B}}_{\widehat{R}}\subset {\bf{H}}^{4}$. It follows from the inner product of the first equation of system (2.1) in ${\bf{H}}$ with $S$, we have that

where ${\bf{H}}_{\frac{1}{2}}$ is the fractional power subspace generated by the sectorial operator $L$. From Theorem 2.2 and Theorem 2.4, we can obtain that $S$ is positive and bounded. Hence, there exists a $C_{1} > 0$ such that

Since ${\bf{H}}_{\frac{1}{2}}\hookrightarrow {\bf{H}}$, so there exists a $C > 0$, such that

As we know that

hence,

By Lemma 2.3, we obtain that

where $\gamma _{1} = 2d_{S}C^{2}$. It follows from the inner product of the second equation of system (2.1) in ${\bf{H}}$ with $V$, we have that

From Theorem 2.2 and Theorem 2.4, we can obtain that $S, V$ are positive and bounded. Hence, there exists a $C_{2} > 0$ such that

Since ${\bf{H}}_{\frac{1}{2}}\hookrightarrow {\bf{H}}$, so there exists a $C > 0$, such that

Hence,

By Lemma 2.3, we obtain that

where $\gamma _{2} = 2d_{V}C^{2}$. It follows from the inner product of the third equation of system (2.1) in ${\bf{H}}$ with $I$, we have that

From Theorem 2.2 and Theorem 2.4, we can obtain that $I, R$ are positive and bounded. Hence, there exists a $C_{3} > 0$ such that

Since ${\bf{H}}_{\frac{1}{2}}\hookrightarrow {\bf{H}}$, so there exists a $C > 0$, such that

Hence,

By Lemma 2.3, we obtain that

where $\gamma _{3} = 2d_{I}C^{2}$. It follows from the inner product of the forth equation of system (2.1) in ${\bf{H}}$ with $R$, we have that

From Theorem 2.2 and Theorem 2.4, we can obtain that there exists a $C_{4} > 0$ such that

Since ${\bf{H}}_{\frac{1}{2}}\hookrightarrow {\bf{H}}$, so there exists a $C > 0$, such that

Hence,

By Lemma 2.3, we obtain that

where $\gamma _{4} = 2d_{R}C^{2}$. Therefore, it implies if $\widehat{R} ^{2} > \max \left\{ \frac{2C_{1}}{\gamma _{1}}, \frac{2C_{2}}{\gamma _{2}}, \frac{2C_{3}}{\gamma _{3}}, \frac{2C_{4}}{\gamma _{4}}\right\}$, then there exists a $t_{0} > 0$, such that for any $t\geq t_{0}$, there has $u\left(t, \varphi \right) \subset {\bf{B}}_{\widehat{R}}$. Thence, ${\bf{B}}_{ \widehat{R}}\subset {\bf{H}}^{4}$ is the absorbing set.

Next we verify Condition (C*). Since $L_{1} = d_{S}\Delta :{\bf{H}} _{1}\rightarrow {\bf{H}}$ is a symmetrical sectorial operator, so the eigenvectors $\left\{ e_{j}\right\} _{j\in \mathbb{N}}$ corresponding to the eigenvalues $\left\{ \lambda _{S, j}\right\} _{j\in \mathbb{N}}$ are the complete orthonormal basis of ${\bf{H}}$, that is for any $S\in {\bf{H}}$,

In addition, $\forall N_{S} > 0, \exists K_{S}\geq 1$ such that $-N_{S}\geq \lambda _{S, j}, \forall j\geq K_{S}+1$.

Let

Then, each $S\in {\bf{H}}$ can be decomposed as

where $P:{\bf{H}}\rightarrow {\bf{H}}_{1}^{K_{S}}$ is the orthogonal projector. There are similar decomposition forms for each $V, I, R\in {\bf{H }}$.

Since $Q\left(t\right)$ has a bounded absorbing set$\ {\bf{B}}_{\widehat{ R}}\subset {\bf{H}}^{4}$, then for any bounded set ${\bf{B}}\subset {\bf{H}}^{4}$, there exists a $t_{\bf{B}} > 0$, without loss of generality, we assume that $t_{\bf{B}} > t_{0}$ such that $\left(S\left(t, S_{0}\left(x\right) \right), V\left(t, V_{0}\left(x\right) \right), I\left(t, I_{0}\left(x\right) \right), R\left(t, R_{0}\left(x\right) \right) \right) \subset {\bf{B}}_{\widehat{R}}$, where

$\forall \left(S_{0}\left(x\right), V_{0}\left(x\right), I_{0}\left(x\right), R_{0}\left(x\right) \right) \in {\bf{B, }}t\geq t_{\bf{B}}$, then

Hence,

It means that $\left\{ \left\Vert PQ\left(t\right) {\bf{B}}\right\Vert _{{\bf{H}}^{4}}\right\} _{t\geq t_{\bf{B}}}$ is bounded.

It follows from the inner product of the first equation of system (2.1) in ${\bf{H}}$ with $S_{2}$, we have that

where

Since $t\geq t_{\bf{B}}$, $\left\Vert S\right\Vert _{\bf{H}}\leq \widehat{R}$. Then

Hence,

By Lemma 2.3,

and $N_{S} > 0$ is arbitrary. Similarly,

and $N_{V}, N_{I}, N_{R} > 0$ are arbitrary. Therefore, we deduce that the Condition (C*) holds. By [53,Theorem 4.1], we obtain that system (2.1) has a global exponential attract set $\mathcal{A}^{\ast }$.

After getting the global exponential attract set, we discuss the stability and uniform persistence of the epidemic disease.

Theorem 3.3. The following statements are valid.

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

in ${\bf{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, V, I, R)$ satisfies

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

Proof. (1) Suppose $\lambda ^{\ast } < 0$. Using the comparison principle, we can show that $I\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 (2.1) that

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

By the comparison principle Lemma 2.1, $I\left(x, t\right) \leq \widetilde{I}\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{I} \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 $I\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 I(x, t)\leq \varepsilon$ for all $x\in \overline{\Omega }, t\geq T$. From the first equation in system (2.1), it is easily observed that $S$ is a super-solution to

and a sub-solution to

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

For system (3.4), we can verify that

It is similar to the proof of Theorem 3.2, we can obtain that system (3.4) exists a global exponential attract set $\mathcal{A}_{w}$. In addition, system (3.4) has a variational structure, the corresponding functional of the variational structure is

where

Then

where $C > 0$, so $L+G$ is a gradient type operator. From [49,Theorem A.2.2], we can prove that

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

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

Thus, our analysis implies that $S\left(\cdot, t\right) \rightarrow S^{0}\left(x\right)$ uniformly on as $t\rightarrow \infty$. Similarly, we can also prove that $V\left(\cdot, t\right) \rightarrow V^{0}\left(x\right)$ uniformly on 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 (2.1). Similar to the proof of conclusion (1), system (3.6) is also a gradient type equation. From [49,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$. Similarly, we can also prove that $V_{-}^{\ast }\left(x\right) > 0$ for all $x\in \Omega$. Next, let us define $\left(I_{-}\left(x, t\right), R_{-}\left(x, t\right) \right) = \left(\varepsilon \varphi ^{\ast }\left(x\right), \varepsilon \psi ^{\ast }\left(x\right) \right)$ where $\varphi ^{\ast }\left(x\right) > > 0, \psi ^{\ast }\left(x\right) > > 0$ are the eigenvalue and eigenvectors in equations (3.2) and $\varepsilon > 0$ is a sufficiently small constant. Substituting $\varepsilon \varphi ^{\ast }, \varepsilon \psi ^{\ast }$ into the third and the forth equations of system (2.1), we know

and

Therefore, $\left(\varepsilon \varphi ^{\ast }, \varepsilon \psi ^{\ast }\right)$ is the sub-solution of the third and the forth equations of system (2.1). We choose $0 < \gamma \left(x\right) < \min \left\{ S_{-}^{\ast }\left(x\right), V_{-}^{\ast }\left(x\right), \varepsilon \varphi ^{\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$, then the positive solution of system (2.1) is globally exponential attractive and the attraction domain is $\mathcal{A}^{\ast }$. From an epidemiological point of view, the disease persists in this situation. The epidemic is globally asymptotically stable or attracted by the global exponential attract set.

4.   Numerical simulation

In this section, we hope to select a specific epidemic disease with some more practical and more definitive data to simulate the stability of our $SVIR$ model with relapse. As we know, HBV epidemic model is a kind of infectious disease model with relapse and vaccination is the most effective way to prevent hepatitis. WHO also releases global hepatitis report every year. From "Global Hepatitis Report, 2017" ^[63], we know that in worldwide, approximately 240 million people have chronic hepatitis B virus infection. Without an expanded and accelerated response, the number of people living with hepatitis B virus is projected to remain at the current, high levels for the next 40-50 years, with a cumulative 20 million deaths occurring between 2015 and 2030. A stepped-up global response can no longer be delayed. From World Health Statistics ^[64], we can find that the data of global natural mortality, the male natural mortality rate is 19% and female natural mortality rate is 12.9% respectively. After a simple calculation we can get the natural mortality rate of the global is 15.95%. In conclusion, we can get the data in Table 2.

According to the data in Table 2, we will first simulate the stability of the constant coefficient model about system (2.1) when $\Lambda, \beta_1, \beta_2, \alpha, \rho, \phi, \mu, \eta _{1}$ and $\eta _{2}$ are assumed to be positive constants. Then we can obtain the stability of disease-free equilibrium of constant coefficient model (Figure 2).

From the numerical simulation, it can be seen that although hepatitis B is still the primary disease endangering human health, however, the global hepatitis B disease has shown a trend of extinction.

In this article, we discuss the dynamics of the epidemic model with a spatial heterogeneous environment, so we then simulate the stability and persistence of the spatial system (2.1). From the system (2.1) we can see that all the parameters are $x-related$ functions, so we choose different functions will directly lead to different stability results. We choose $\beta _{1}\left(x\right) = \frac{1}{1+x}, \mu \left(x\right) = e^{-2x}, \phi \left(x\right) = e^{-x}, \rho \left(x\right) = 0.4\sin x, \eta _{1}\left(x\right) = 0.2, \eta _{2}\left(x\right) = 0.1$ and select other parameters from Table 2. We can show the global stability of non-constant disease-free equilibrium $E^{0}\left(x\right)$ of system (2.1) (Figure 3).

We choose $\mu \left(x\right) = e^{-x}, \phi \left(x\right) = \frac{1}{1+x}, \beta _{1}\left(x\right) = 0.62, \beta _{2}\left(x\right) = 0.33e^{-x}$ and select other parameters from Table 2. We can obtain that the spatial heterogeneity epidemic is global stability (Figure 4).

We choose $\mu \left(x\right) = e^{-x}, \phi \left(x\right) = \frac{1}{1+x}, \rho \left(x\right) = 0.4\sin x, \beta _{1}\left(x\right) = 0.62, \beta _{2}\left(x\right) = 0.33e^{-x}$ and select other parameters from Table 2. We can clearly see that the spatial heterogeneity epidemic is persists uniformly (Figure 5).

In Figure 5, we can find that the image of the solution is turbulent, and the amplitude decreases with time. Finally, the image of the turmoil is controlled in range. And this range is the scope of the global exponential attract set. From the comparison of Figure 4 and Figure 5, we can find that spatial heterogeneity has a great impact on the spread of disease, and different parameters lead to changes in long-term dynamic behavior.

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

5.   Conclusions and discussion

We have formulated a reaction-diffusion $SVIR$ model with relapse in the spatially heterogeneous environment. We first prove the positivity of solution by using the comparison principle of monotone nonlinear reaction-diffusion dynamical system and comparison principle of parabolic equation. Most of the models in the previous work are non-coupled models or systems that contain only one diffusion equation. The usual way to study these two problems is to study one of the equations in the model first (such as eigenvalue problems or the existence of solutions), and then to get the properties of several other equations. However, the above methods do not applicable to the coupling system of the four diffusion equations we have established in this paper. In addition, we prove that the disease-free equilibrium is globally asymptotically stable and the endemic is persisting uniformly by the global exponential attract set theory. From our proof, we can find that $\lambda ^{\ast }$ has threshold characteristics. From the numerical simulation we can find the effect of spatial heterogeneity on the diffusion of disease is very large. Simulation results show that the best way to eliminate the disease is to control the contact rate. Reduce the relapse rate or increase the vaccination rate can make the disease tends to a stable state. Choosing different spatial coefficients can produce a fundamental change of the global dynamics of the system.

Acknowledgements

The authors would like to thank the anonymous referees for their useful and valuable suggestions. The first author was supported by the Fundamental Research Funds for the Central Universities (JUSRP11949), NSF of China (11671180, 11671181) and the third author was supported by the science and technology foundation of Suqian(Z201444), China.

Conflict of interest

The authors declare there is no conflict of interest.