Global behavior of a multi-group SEIR epidemic model with age structure and spatial diffusion

1.   Introduction

Since the pioneering work of Kermack and McKendrick ^[1], many mathematical models have been proposed attempting to gain a better understanding of disease transmission, especially for the control strategy and dynamical behavior of infectious diseases ^{[2,3,4,5,6,7,8]}. Simple models with assumption that individuals are well mixed, which implies each individual has the same probability to be infected, are beneficial in that one can obtain analytical results easily but may be lack of realism. Epidemic models with population structures, like age, sex and patch (such as communities, cities, or counties), may be a more realistic way to describe complex disease dynamics. As a matter of fact, the total population should be classified into different groups and the vital epidemic parameters should vary among different population groups. In addition, at different age stages, the effects of infectious transmission are various, which is another important and key factor that needs necessarily to be included in modeling this infectious transmission process. Thus, considering multi-group and age structure in epidemic models is very necessary and reasonable. Some recent developments on the transmission dynamics of multi-group and age structured epidemic models have been discussed in ^{[8,9,10,11,12,13]}.

Since that the population distribute heterogeneously in different spatial location in the real life and they will move or diffuse for many reasons, in epidemiology, there is increasing evidence that environmental heterogeneity and individual motility have significant impact on the spread of infectious diseases ^[14,15]. In recent years, global behavior of spatial diffusion systems, which are suitable for diseases such as the rabies and the Black Death, has been attracted extensive attention of researchers and has been one of the hot topics ^{[16,17,18,19,20,21,22,23,24,25]}. Among these works, few take age structure or multi-group into consideration. Yang et al. ^[24] proposed a novel model incorporated with both age-since-infection and spacial diffusion of brucellosis infection, and the basic reproduction number and global behaviors of this system were completely investigated. Fitzgibbon et al. ^[19] considered a diffusive epidemic model with age structure where the disease spreads between vector and host populations. Then, the existence of solutions of the model was studied based on semigroup theory and the asymptotic behavior of the solution was analyzed. Luo et al. ^[21] incorporated spatial heterogeneity in $n$-group reaction-diffusion SIR model with nonlinear incidence rate to investigate the global dynamics of the disease-free and endemic steady states for this model. Zhao et al. ^[25] modeled host heterogeneity by introducing multi-group structure in a time delay SIR epidemic model and showed that basic reproduction number determines the existence of traveling waves of this system. To determine how age structure, multi-group population and diffusion of individuals affect the consequences of epidemiological processes, Ducrot et al. ^[18] formulated a multi-group age-structured epidemic model with the classical Fickian diffusion and studied the existence of travelling wave solutions for this model.

To the best of our knowledge, epidemic models established by researchers except for ^[18] only include one or two characteristics of multi-group, age structure and spatial diffusion. All the three characteristics are incorporated into epidemic model in ^[18], however, this model does not include the class of latent individuals. For some epidemic diseases like malaria, HIV/AIDS and West Nile virus, latent individuals may take days, months, or even years to become infectious. Moreover, the travel of latent individuals showing no symptoms can spread the disease geographically which makes disease harder to control. Motivated by the above discussion, in this paper, we investigate a diffusive version of multi-group epidemic system with age structure which is generalization of the model studied in ^[26] for the first time to allow for individuals moving around on the spatial habitat $x\in\Omega\subset\mathbb{R}^n$ with smooth boundary $\partial\Omega$. The organization of this paper is as follows. Firstly, we present our model in the next section. In section 3, some preliminaries including the positivity, boundedness, existence and uniqueness of solution, and the existence of compact global attractor of the associated solution semiflow, are established. In section 4, the sufficient conditions on the existence and global stability of disease-free and endemic steady states are stated and proved. In section 5, we conduct numerical simulations to illustrate the validity of our theoretical results. In section 6, a brief conclusion is given.

2.   The model

In 2015, Liu et al. ^[26] introduced age-of-latent and age-of-relapse into epidemic model which is appropriate for diseases such as tuberculosis and herpes virus infection. For these diseases, latent individuals may take days, months, or even years to become infectious and the treatment efficacy may decline with time for recovered individuals then cause recurrence of disease. In order to study the global dynamics o for these diseases, they formulated the following SEIR epidemic system with continuous age-dependent latency and relapse

for $t\geq0$ and with initial conditions

for $a, b\geq0$, where $L_+^1(0, +\infty)$ is the space of functions on $(0, +\infty)$ that are nonnegative and Lebesgue integrable. At time $t$, the densities of susceptible individuals, latent individuals with latent age $a$, infectious individuals and removed individuals with relapse age $b$ are denoted by $S(t)$, $e(t, a)$, $I(t)$, $r(t, b)$, respectively. $\sigma(a)$ and $\gamma(b)$ denote the conversional rate from the latent class and the relapse rate in the removed class, which depend on age $a$ and age $b$, respectively. Furthermore, $\beta$ is the transmission rate of the disease between susceptible and infectious individuals, $\Lambda$ is the density of the recruitment into the susceptible class (including the births and immigration), $\mu$ is the natural death rate of all individuals, $\delta_1$ and $\delta_2$ are the additional death rate induced by the infectious diseases, and $c$ is the recovery rate from the infectious class. All parameters are assumed to be positive.

It is clear that the variations of different epidemic parameters between or within different groups can be well realized according to the description of multi-group epidemic models. Hence, Liu and Feng ^[27] extended model (2.1) to the situation in which the population is divided into $n$ groups according to different contact patterns and derived the following multi-group SEIR epidemic model

for $t\geq0$ and with initial conditions

for $a, b > 0$. $\Lambda_k$, $\mu_k$ and $c_k$ denote the recruitment rate of the susceptible class, the per-capita natural death rate and the recovery rate from the infectious class in group $k$, respectively. $\beta_{kj}$ denotes the transmission rate of the disease between susceptible individuals in group $k$ and infectious individuals in group $j$. $\delta_{1k}$ and $\delta_{2k}$ denote the additional death rates of exposed and infectious individuals induced by the infectious diseases in group $k$, respectively. $\sigma_k(a)$ denotes the conversional rate from the latent class in group $k$, which depends on age $a$ and $\gamma_k(b)$ denotes the relapse rate from the removed class into the infectious class in group $k$, which depends on age $b$.

Spatial diffusion is an intrinsic characteristic for investigating the roles of spatial heterogeneity on diseases mechanisms and transmission routes and can lead to rich dynamics. Based on this fact, we generalize (2.2) by taking account of the case that individuals move or diffuse around on the spatial habitat $x\in\Omega\subset\mathbb{R}^n$ with smooth boundary $\partial\Omega$. Let $S_k(t, x)$ and $I_k(t, x)$ be the densities of susceptible individuals and infectious individuals at time $t$ and location $x\in\Omega$ in group $k$, respectively, where the habitat $\Omega$ is bounded and connected. And let $e_k(t, a, x)$ and $r_k(t, b, x)$ denote the densities of individuals in the latent class with age $a$ and the removed class with age $b$ at time $t$ and location $x$ in group $k$, respectively. Hence, the $n$-group diffusive SEIR epidemic model with age-dependent latent and relapse has the following form

for $x\in\Omega$, $a, b\in \mathbb{R}_+ = (0, +\infty)$, with the homogeneous Neumann boundary conditions

and initial functions

$d_{1k}$, $d_{2k}$, $d_{3k}$, $d_{4k}$ denote the diffusion coefficients of susceptible individuals, exposed individuals, infectious individuals and removed individuals in group $k$, respectively. And the other parameters have the same biological meanings as in (2.2). The homogeneous Neumann boundary conditions imply that there is no population flux across the boundary $\partial\Omega$.

We define the functional spaces $X = C(\bar{\Omega}, \mathbb{R})$ and $Y = L^1(\mathbb{R}_+, X)$ for model (2.3) equipped, respectively, with the norms

for $\phi\in X$, $\varphi\in Y$. The positive cones are denoted by $X_+$ and $Y_+$. In addition, we define a vector space $Z = (C([0, T], X))^{2n}$ with the norm

Throughout this paper, for convenience, we always denote $S = (S_1, S_2, ..., S_n)$, $e = (e_1, e_2, ..., e_n)$, $I = (I_1, I_2, ..., I_n)$, $r = (r_1, r_2, ..., r_n)$, and $S^0 = (S^0_1, S^0_2, ..., S^0_n)$, $e^0 = (e^0_1, e^0_2, ..., e^0_n)$, $I^0 = (I^0_1, I^0_2, ..., I^0_n)$, $r^0 = (r^0_1, r^0_2, ..., r^0_n)$. We also denote $(y_1, y_2, ..., y_n)^T > (z_1, z_2, ..., z_n)^T$ as $y_i > z_i$ for all $i = 1, 2, ..., n$. For each $i = 1, 2, 3, 4$, we suppose that $T_{ik}:C(\bar{\Omega}, \mathbb{R})\rightarrow C(\bar{\Omega}, \mathbb{R})$ is the $C_0$ semigroup generated by $d_{ik}\Delta$ subjects to the Neumann boundary condition in group $k$. From subsection 2.1 in ^[28], we have

for all $t > 0$ and $\phi\in C(\bar{\Omega}, \mathbb{R})$, where $\Gamma_{ik}(t, x, y)$ is the Green function. We have that $T_{ik}$, $i = 1, 2, 3, 4$, $k = 1, 2, ..., n$ are compact and strongly positive for each $t > 0$ by the Corollary 7.2.3 in ^[29]. Integrating the second equation in model (2.3) along the characteristic line $t-a = c$, where $c$ is a constant, we obtain

where $\pi_{1k}(a) = e^{-\int_0^a[\mu_k+\delta_{1k}+\sigma_k(s)]ds}$. Similarly,

where $\pi_{2k}(b) = e^{-\int_0^b[\mu_k+\gamma_k(s)]ds}$. To study the asymptotic behaviors of the dynamics of model (2.3), we require the following assumptions on the model parameters.

Assumption 2.1. For each $k, j = 1, 2, ..., n$,

$(H_1) \; d_{1k}, d_{2k}, d_{3k}, d_{4k}, \Lambda_k, \mu_k, \delta_{1k}, \delta_{2k}, c_k > 0$.

$(H_2) \; \beta_{kj}\in\mathbb{R}_+$, and the $n$-dimensional square matrix $(\beta_{kj})_{n\times n}$ is irreducible.

$(H_3) \; \sigma_k(\cdot), \gamma_k(\cdot)\in L^{\infty}(\mathbb{R}_+, \mathbb{R}_+)$, $\bar\sigma_k > 0$ and $\bar\gamma_k > 0$, where $\bar{\sigma}_k: = \mathop{\rm{ess\; sup}}\limits_{a\in\mathbb{R}_+}\sigma_k(a)$, $\bar{\gamma}_k: = \mathop{\rm{ess\; sup}}\limits_{b\in\mathbb{R}_+}\gamma_k(b).$

3.   Preliminaries

We define $A_k(t, x) = e_k(t, 0, x)$, $B_k(t, x) = r_k(t, 0, x)$ for $(t, x)\in\mathbb{R}_+\times\Omega$, and let $A = (A_1, A_2, ..., A_n)$, $B = (B_1, B_2, ..., B_n)$ and $C = (A, B)$. Positiveness of the solutions of model (2.3) is given below.

Theorem 3.1. Suppose that there exists a solution $(S(t, \cdot), e(t, \cdot, \cdot), I(t, \cdot), r(t, \cdot, \cdot))\in X^n\times Y^n\times X^n\times Y^n$ of (2.3) corresponding to $(S^0, e^0, I^0, r^0)\in X_+^n\times Y_+^n\times X_+^n\times Y_+^n$ with an interval of existence $[0, T]$, $T > 0$. Then

for all $t\in[0, T]$.

Proof. From the first equation of (2.3), we have

where $\hat{F}_{S_k}(t, x) = e^{-\int_{0}^t\big[\mu_k+{\sum\limits_{j = 1}^n}{\beta_{kj}I_j(\tau, x)}\big]d\tau}\int_\Omega\Gamma_{1k}(t, x, y)S_k^0(y)dy$. The positivity of $\Lambda_k$ and $S_k^0$ ensures $S_k(t, x) > 0$ for each $(t, x)\in[0, T]\times\Omega$. The positivity of $C$ which means the positivity for $A_k$ and $B_k$, $k = 1, 2, ..., n$ is established by constructing Picard sequences as follows.

Solving equation $I_k$ for system (2.3), we have

where $F_{I_k}(t, x) = e^{-(\mu_k+\delta_{2k}+c_k)t}\int_\Omega\Gamma_{3k}(t, x, y)I_k^0(y)dy.$ For $(t, x)\in[0, T]\times\Omega$, by (2.4) and (2.5), we obtain

where $F_{A_k}(t, y) = \int_0^\infty\sigma_k(b+t)\frac{\pi_{1k}(b+t)}{\pi_{1k}(b)}\int_\Omega\Gamma_{2k}(b+t, y, z)e_k^0(b, z)dzdb$, and

where $F_{B_k}(t, y) = \int_0^\infty\gamma_k(b+t)\frac{\pi_{2k}(b+t)}{\pi_{2k}(b)}\int_\Omega\Gamma_{4k}(b+t, y, z)r_k^0(b, z)dzdb$. From (3.1)–(3.3) and the definitions of $A_k$ and $B_k$, we have

and

Let

and

Then it is obvious that $A_k^{(0)}(t, x) > 0$, $B_k^{(0)}(t, x) > 0$. Now we assume that $A_k^{(m)}(t, x) > 0$, $B_k^{(m)}(t, x) > 0$ ($m\in\mathbb{N}$) for $e^0_k > 0$, $r^0_k > 0$ and $(t, x)\in[0, T]\times\Omega$. Then

and

From the positivity of $\beta_{kj}$, $\sigma_k$ and $\gamma_k$, together with the positivity of $\Gamma_{2k}$, $\Gamma_{3k}$ and $\Gamma_{4k}$, it follows that

and

which lead to $C^{(1)}(t, x)-C^{(0)}(t, x) > 0$ for $(t, x)\in[0, T]\times\Omega$. We assume that $C^{(m)}(t, x)-C^{(m-1)}(t, x) > 0$ for all $m\geq2$, that is, $A_k^{(m)}(t, x)-A_k^{(m-1)}(t, x) > 0$ and $B_k^{(m)}(t, x)-B_k^{(m-1)}(t, x) > 0$, $k = 1, 2, ..., n$. Then,

and

Hence applying mathematical induction, we show that the sequence $\{C^{(m)}\}_0^\infty$ is monotonically increasing.

Next, applying the contraction mapping principle, we show the sequence $\{C^{(m)}\}_0^\infty$ converges to $C(t, x)$ for any $(t, x)\in[0, T]\times\Omega$ as $m$ approaches infinity. To this end, we define a variable

By the definitions of $A_k^{(m)}$ and $B_k^{(m)}$, we have

and

For any $m\in\mathbb{N}$,

where $\hat{S}_k = \max\limits_{t\in[0, T]}|S_k(t, \cdot)|_X$, and

Hence,

where $K_\lambda = \max\left\{M_\lambda, N_\lambda\right\}$, $M_\lambda = \max\limits_k\bigg\{{\sum\limits_{j = 1}^n}\frac{{\beta_{kj}}\hat{S}_k(\bar{\sigma}_j+\bar{\gamma}_j)}{\lambda^2}\bigg\}$, $N_\lambda = \max\limits_k\left\{\frac{c_k(\bar{\sigma}_k+\bar{\gamma}_k)}{\lambda^2}\right\}$. Therefore, for any $m_1 > m_2$, $m_1, m_2\in\mathbb{N}$,

We choose $\lambda$ sufficiently large such that ${\sum\limits_{j = 1}^n}\frac{{\beta_{kj}}\hat{S}_k(\bar{\sigma}_j+\bar{\gamma}_j)}{\lambda^2} < 1$ and $\frac{c_k(\bar{\sigma}_k+\bar{\gamma}_k)}{\lambda^2} < 1$ for all $k = 1, 2, ..., n$, then $K_\lambda < 1$. Hence, $\|\hat{C}^{(m_1)}-\hat{C}^{(m_2)}\|_Z\rightarrow0$ as $m_2\rightarrow\infty$ which implies that $\hat{C}^{(m)}\rightarrow\hat{C}$ and thus $C^{(m)}\rightarrow C$ as $m\rightarrow\infty$. Furthermore, we have $A_k^{(m)}\rightarrow A_k$ and $B_k^{(m)}\rightarrow B_k$ for $k = 1, 2, ..., n$ as $m\rightarrow\infty$. Since sequence $\{C^{(m)}\}_0^\infty$ is monotonically increasing, we obtain $A_k$ and $B_k$ are positive for $k = 1, 2, ..., n$.

By (2.4) and (2.5), together with the positivity of $e_k^0$, $r_k^0$, $A_k$ and $B_k$, we conclude that $e_k(t, a, x)$ and $r_k(t, a, x)$ are positive. For the positivity of $I_k$, we prove this by contradiction. Suppose that there exist $x_0\in\Omega$ and $t_0 = \inf\{t\in\mathbb{R}_+|I_k(t, x_0) = 0\}$ such that

By the third equation of system (2.3), we can easily obtain

This leads to a contradiction. Hence, for any $t\in[0, T]$, we have $(S(t, \cdot), e(t, \cdot, \cdot), I(t, \cdot), r(t, \cdot, \cdot))\in X_+^n\times Y_+^n\times X_+^n\times Y_+^n$.

Let $N_k(t) = \int_\Omega\left[S_k(t, x)+\int_0^\infty e_k(t, a, x)da+I_k(t, x)+\int_0^\infty r_k(t, b, x)db\right]dx$ denotes the total population at time $t$ in group $k$ and region $\Omega$.

Theorem 3.2. If

for all $t > 0$, $x\in\Omega$, the region $\Pi$ defined by

is positively invariant for system (2.3).

Proof. Following condition (3.4) and the equations of system (2.3), we have

Noting the Neumann boundary conditions of system (2.3) and using the Gauss formula, we derive

It follows that

Thus if $N_k(t) > \frac{\Lambda_k}{\mu_k}|\Omega|$, then $\frac{dN_k(t)}{dt} < 0$. Moreover, we observe the ordinary differential equation

with general solution

where $N_k(0)$ means the initial value of total population in group $k$ and region $\Omega$. By applying the standard comparison theorem, we have for all $t\geq0$,

Hence, $\Pi$ is positive invariant for system (2.3).

The existence and uniqueness of the solution of model (2.3) follow from Banach-Picard fixed point theorem.

Theorem 3.3. Let initial functions satisfy $(S^0, e^0, I^0, r^0)\in X_+^n\times Y_+^n\times X_+^n\times Y_+^n$. Then the system (2.3) has a unique solution $(S(t, \cdot), e(t, \cdot, \cdot), I(t, \cdot), r(t, \cdot, \cdot))\in X_+^n\times Y_+^n\times X_+^n\times Y_+^n$ for $t\in[0, T]$.

Proof. Solving equation $S_k$ for system (2.3), we have

for $(t, x)\in[0, T]\times\Omega$, where $F_{S_k}(t, x) = e^{-\mu_kt}\int_\Omega\Gamma_{1k}(t, x, y)S_k^0(y)dy.$ From (3.1)–(3.5) and the definitions of $A_k$ and $B_k$, we have

and

where $\mathcal{F}_{1k}$, $\mathcal{F}_{2k}:Z\rightarrow C([0, T], X)$ are nonlinear operators for each $k = 1, 2, ..., n$. For the sake of convenience, we define for each $(t, x)\in[0, T]\times\Omega$,

Then

For any $C$, $\bar{C}\in Z$, we set $\tilde{C} = C-\bar{C}$. Then, from the positivity of $A_k$ and $B_k$ proved in Theorem 3.1, we have

and

where

Denote

and

Clearly, we can choose $T$ small enough such that $m_{1k}(T) < 1$ and $m_{2k}(T) < 1$ for all $k = 1, 2, ..., n$. Consequently, we have $m(T) < 1$. Then

Hence, applying contraction operator theorem ^[30], we conclude that $\mathcal{F}$ has a unique fixed point $C = (A_1, A_2, ..., A_n, B_1, B_2, ..., B_n)$. From Theorem 3.1, (2.4) and (2.5), together with $A_k(t, x) = {\sum\limits_{j = 1}^n}{\beta_{kj}S_k(t, x)I_j(t, x)}$ and $B_k(t, x) = c_kI_k(t, x)$, we derive the existence and uniqueness of the solution $(S(t, \cdot), e(t, \cdot, \cdot), I(t, \cdot), r(t, \cdot, \cdot))\in X_+^n\times Y_+^n\times X_+^n\times Y_+^n$ for system (2.3).

To further establish the global existence of the solution of system (2.3), we need the following lemma.

Lemma 3.1. ^[31]. The following problem

admits a unique positive steady state $\omega^* = \frac{\Lambda}{\mu}$, which is globally attractive in $X$.

Theorem 3.4. Let initial functions satisfy $(S^0, e^0, I^0, r^0)\in X_+^n\times Y_+^n\times X_+^n\times Y_+^n$. Then the system (2.3) has a unique solution $(S(t, \cdot), e(t, \cdot, \cdot), I(t, \cdot), r(t, \cdot, \cdot))\in X_+^n\times Y_+^n\times X_+^n\times Y_+^n$ for $t\in\mathbb{R}_+$.

Proof. To extend the domain of existence from $t\in[0, T]$ to $t\in\mathbb{R}_+$, it suffices to show that the solution does not blow up in finite time. In fact, by Theorem 3.1, we have that

for all $t > 0$ and $x\in\Omega$. From Lemma 3.1 and the comparison principle, we get that $S_k(t, x)$ is bounded above by the upper solution $\frac{\Lambda_k}{\mu_k}$.

We now claim that $e_k(t, a, x) < +\infty$ for all $t > 0$, $a > 0$, $x\in\Omega$ and $k = 1, 2, ..., n$. From (2.4), it is sufficient to show that $e_k(t, 0, x) < +\infty$ for all $t > 0$ and $x\in\Omega$. Suppose on the contrary that there exist $t_e > 0$ and $x_e\in\Omega$ such that

We then have from (3.5) that

which implies that $S_k(t, x_e)$ is negative in the neighborhood of $t_e$. This contradicts to the positivity of $S_k$, which has been proved in Theorem 3.1. Furthermore, from $e_k(t, 0, x) = {\sum\limits_{j = 1}^n}{\beta_{kj}S_k(t, x)I_j(t, x)}$ and $r_k(t, 0, x) = c_kI_k(t, x)$ in (2.4), we obtain $I_k(t, x) < +\infty$ and $r_k(t, 0, x) < +\infty$. And from (2.5), we get $r_k(t, a, x) < +\infty$ for all $t > 0$, $a > 0$, $x\in\Omega$ and $k = 1, 2, ..., n$. Thus, blow up never occurs, and we obtain a solution $(S(t, \cdot), e(t, \cdot, \cdot), I(t, \cdot), r(t, \cdot, \cdot))\in X_+^n\times Y_+^n\times X_+^n\times Y_+^n$ for $t\in\mathbb{R}_+$.

Theorem 3.5. The solution semiflow $\Phi(t) = \left(S(t, \cdot), e(t, \cdot, \cdot), I(t, \cdot), r(t, \cdot, \cdot)\right):X^n_+\times Y^n_+\times X^n_+\times Y^n_+\backslash(\bf{0}, \bf{0}, \bf{0}, \bf{0})\rightarrow X^n_+\times Y^n_+\times X^n_+\times Y^n_+\backslash(\bf{0}, \bf{0}, \bf{0}, \bf{0})$ of model (2.3) has a compact and global attractor as condition (3.4) holds.

Proof. According to Theorem 3.2, we know that system (2.3) is ultimately bounded, which implies that solution semiflow $\Phi(t)$ is point dissipative on $X^n_+\times Y^n_+\times X^n_+\times Y^n_+\backslash(\bf{0}, \bf{0}, \bf{0}, \bf{0})$. By Theorem 2.6 in ^[15], we can get that $\Phi(t)$ is compact for any $t > 0$. Thus, from Theorem 3.4.8 in ^[32], we further obtain that $\Phi(t)$ has a compact and global attractor in $X^n_+\times Y^n_+\times X^n_+\times Y^n_+\backslash(\bf{0}, \bf{0}, \bf{0}, \bf{0})$.

4.   Global stability

4.1.   Existence of equilibria

It is easy to see that model (2.3) has a unique disease-free steady state $\bar{E} = (\bar{S}, \bf{0}, \bf{0}, \bf{0})^T$ where $\bar{S} = (\bar{S}_1, \bar{S}_2, ..., \bar{S}_n)^T$ and $\bar{S}_k = \frac{\Lambda_k}{\mu_k}$, $k = 1, 2, ..., n.$

We denote $E^* = (S^*, e^*(\cdot), I^*, r^*(\cdot))^T$ as the space-independent endemic steady state of (2.3), where $S^* = (S_1^*, S_2^*, ..., S_n^*)^T$, $e^*(\cdot) = (e_1^*(\cdot), e_2^*(\cdot), ..., e_n^*(\cdot))^T$, $I^* = (I_1^*, I_2^*, ..., I_n^*)^T$ and $r^*(\cdot) = (r_1^*(\cdot), r_2^*(\cdot), ..., r_n^*(\cdot))^T$. Then, $E^*$ satisfies

We denote $f_k(I^*) = {\sum\limits_{j = 1}^n}{\beta_{kj}I_j^*}:\mathbb{R}_+^n\rightarrow\mathbb{R}_+$, $h_k(I^*) = \frac{\Lambda_k}{\mu_k+f_k(I^*)}:\mathbb{R}_+^n\rightarrow\mathbb{R}_+$, $k = 1, 2, ..., n$ for brevity. By solving the Eq (4.1), we get

where $\pi_{1k}(\cdot)$ and $\pi_{2k}(\cdot)$ are given in (2.4) and (2.5), respectively. Substituting (4.2) into the third equation in (4.1) and rearranging it, we have

where $L_k = {\int_0^\infty}{\sigma_k(a)\pi_{1k}(a)da}$ and $P_k = {\int_0^\infty}{\gamma_k(b)\pi_{2k}(b)db}$. Let us define a matrix-valued function $M(x)$ on $\mathbb{R}^n$ to $\mathbb{R}^{n\times n}$, where $M(x)_{ij} = \frac{L_i\beta_{ij}h_i(x)}{(\mu_i+\delta_{2i}+c_i)-c_iP_i}$, $x = (x_1, x_2, ..., x_n)^T\in\mathbb{R}^n$. Then, (4.3) is equivalent to

On the existence of the endemic equilibrium $E^*$ of system (2.3), we prove the following theorem.

Theorem 4.1. Let $M^0 = \left(\frac{{\bar{S}_iL_i\beta_{ij}}}{(\mu_i+\delta_{2i}+c_i)-c_iP_i}\right)_{n\times n}$. If $\rho(M^0) > 1$, where $\rho(M^0)$ represents the spectral radius of $M^0$, then system (2.3) has a space-independent steady state $E^*$.

Proof. From (4.4), we only need to show that the nonlinear operator $\tilde{M}(x): = M(x)x$, $x\in\mathbb{R}_+^n$, has at least one positive fixed point $x^*\in\mathbb{R}_+^n$. We define $\|\tilde{M}(x)\|: = \max\limits_{1\leq k\leq n}|\tilde{M}(x)_k|$, where $\tilde{M}(x)_k$ denotes the $k$-th entry of vector $\tilde{M}(x)$. Then, $\tilde{M}(x)$ is monotone increasing with respect to $x\in\mathbb{R}_+^n$ and uniformly bounded above by $\max\limits_{1\leq k\leq n}\left|\frac{{\Lambda_kL_k}}{(\mu_k+\delta_{2k}+c_k)-c_kP_k}\right|$.

It is obvious that $\tilde{M}(\bf{0}) = \bf{0}$ and ${M}(\bf{0})$ is the strong Fréchet derivative of $\tilde{M}(\cdot)$ at the origin. Since $M(\bf{0}) = M^0$, we have $\rho(M(\bf{0})) > 1$. Thus, it follows from the Perron-Frobenius theorem (see ^[33]) that $\rho(M(\bf{0}))$ is a simple eigenvalue of $M(\bf{0})$ corresponding to a strictly positive eigenvector, and there exists no nonnegative eigenvector of $M(\bf{0})$ corresponding to eigenvalue 1. Hence, we apply Theorem 4.11 of ^[34], to conclude that operator $\tilde{M}(\cdot)$ has at least one positive fixed point $x^*\in\mathbb{R}_+^n$.

4.2.   Global stability of equilibria

Theorem 4.2. For all $k = 1, 2, ..., n$, if

where $\alpha_k = \frac{{\sum\limits_{j = 1}^n}\bar{S}_j\beta_{jk}}{\mu_k+\delta_{2k}}$, $\pi_{1k}(\cdot)$ and $\pi_{2k}(\cdot)$ are given in (2.4) and (2.5), respectively, then the disease-free steady state $\bar{E}$ is globally asymptotically stable.

Proof. We construct a Lyapunov function

where

and

Taking the derivation of $V_{1k}(t, x)$ along the trajectory of (2.3) with respect to $t$, we have

Recalling (2.4), we can rewrite $V_{2k}(t, x)$ as follows

Thus, we calculate $\frac{\partial V_{2k}}{\partial t}$ along the solution of system (2.3) and get

Similarly,

From (2.5), we get

Thus,

Hence, combining (4.7)–(4.10), we calculate the derivative of $V_{k}(t)$ along the solution trajectory of (2.3) as

Using (4.6), we yield that

Thus, from (4.5), we have

which implies the global asymptotic stability of disease-free steady state $\bar{E}$.

Theorem 4.3. If $\rho(M^0) > 1$, $\lim\limits_{a\rightarrow+\infty}e_k(t, a, x) = 0$ and $\lim\limits_{b\rightarrow+\infty}r_k(t, b, x) = 0$ for all $t > 0$, $x\in\Omega$, then the space-independent steady state $E^*$ is globally asymptotically stable.

Proof. We define $g(x) = x-1-\ln x$, clearly, $g(x)$ is always positive for $x > 0$ and $g'(x) = 1-\frac{1}{x}$. Consider a matrix $D = (\bar{\beta}_{kj})_{n\times n}$ with entry $\bar{\beta}_{kj} = \beta_{kj}L_kS_k^*I_j^*$ and a digraph $\mathcal{G} = (U, H)$ which contains a set $U = \{1, 2, ..., n\}$ of vertices and a set $H$ of arcs $(k, j)$ leading from initial vertex $k$ to terminal vertex $j$, then, we denote a weighted digraph as $(\mathcal{G}, D)$ for which each arc $(j, k)$ is assigned a positive weight $\bar{\beta}_{kj}$. Furthermore, we denote $\bar{D}$ as the Laplacian matrix of matrix $(\mathcal{G}, D)$. Then, the irreducibility of matrix $(\beta_{kj})_{n\times n}$ implies that $\bar{D}$ is also irreducible. Let $q_k$ denote the cofactor of $k$-th diagonal element of $\bar{D}$. And we construct a Lyapunov function as the following form

where

and

For convenience, we denote $J_k(t, x) = {\sum\limits_{j = 1}^n}\beta_{kj}I_j(t, x)$ and $J_k^* = {\sum\limits_{j = 1}^n}\beta_{kj}I_j^*$. Taking the derivation of $W_{1k}$ along the trajectory of (2.3) with respect to $t$, we have

Calculating the derivative of $W_{2k}$ along the solution of system (2.3) yields

Similarly,

and

Hence, combining (4.12)–(4.15), we calculate the derivative of $W_{k}(t)$ along the solution trajectory of (2.3) as

By Theorem 2.3 in ^[35], the following identity holds

Hence, together with the property of $g$, we have $\frac{dW}{dt}\leq0$. Furthermore, the strict equality holds only if $S_k(t, x) = S_k^*$, $e_k(t, a, x) = e_k^*(a)$, $I_k(t, x) = I_k^*$, and $r_k(t, b, x) = r_k^*(b)$. Thus, $T^* = \left\{E^*\right\}\subset\Omega$ is the largest invariant subset of $\left\{(S, e, I, r):\frac{dW}{dt} = 0\right\}$, and the Lyapunov-LaSalle invariance principle implies that the endemic equilibrium $E^*$ is globally asymptotically stable when $\rho(M^0) > 1$.

5.   Numerical simulations

In this section, we present numerical examples to demonstrate the validity and applicability of our main results, Theorems 4.2 and 4.3. For simplicity, we consider the case of two groups, a normalized 1-dimensional space ($\Omega$ = [0, 3]) and a normalized maximum age ($a$ = 2). Firstly, we fix some parameter values as follows,

Then, we take the initial functions as

Example 5.1. If we choose other parameters as

then we have

In this case, from Theorem 4.2, we expect the disease-free steady state $\bar{E}$ to be globally asymptotically stable. In fact, in Figure 1, the infective population $I_1(t, x)$ and $I_2(t, x)$ converge to zero over time.

Example 5.2. If we choose other parameters as

then we have $\rho(M^0)\approx3.89217 > 1$. In this case, from Theorem 4.3, we expect the space-independent steady state $E^*$ to be globally asymptotically stable. In fact, in Figure 2, the infective population $I_1(t, x)$ and $I_2(t, x)$ converge to the positive distribution over time.

6.   Conclusions

In this paper, as an additional structure of the system, we focus on the spatial diffusion of the population. Models with spatial diffusion allow individuals to move to adjacent positions through a random walk process, this is a key factor in considering the geographical spread of infectious diseases. Firstly, we propose the $n$-group diffusive SEIR epidemic model with age-dependent latent and relapse, it is a generalization of the model in ^[27] to a spatially diffusive system. Then, we investigate the positivity, boundedness, existence and uniqueness of solution and the existence of compact global attractor of the associated solution semiflow for this system. For these results, we use the method of constructing Picard sequences, Banach-Picard fixed point theorem and theories of partial functional differential equations. Thereafter, we establish the existence of disease-free and endemic steady states based on Perron-Frobenius theorem. we utilize appropriate Lyapunov functionals, graph-theoretical results and the LaSalle's invariance principle to prove the global stability of disease-free and endemic steady states. Thus, we presented the results of numerical simulations to verify the validity of our main theorems. This is important from the viewpoint of applications.

In this epidemic model, we are concerned with two kinds of spatial heterogeneity: the patch structure and spatial diffusion. Furthermore, age-of-latent and age-of-relapse are included into the epidemic model which is appropriate for diseases such as tuberculosis and herpes virus infection. Dynamical results obtained in this paper provide theoretical foundation for seeking effective measures to prevent, control and study disease transmission.

The expressions of basic reproduction number and endemic steady state depends on space are not analyzed in this paper owing to the complexity of model. In addition, how to improve the sufficient conditions that ensure the stabilities of steady states and make them be depended on basic reproduction number is also need to investigate. We leave these issues for future research.

Acknowledgments

The authors are very grateful to the editor and the anonymous referees for their valuable comments and suggestions, which helped us to improve the presentation of this work significantly.

This work was supported by the National Natural Science Foundation of China (Grant Number: 11971329).

Conflict of interest

The authors declare there is no conflict of interest.