**2. An Age-Structured SEIRS Model with Time Delay**

Motivated by the referred works [28], we discuss two age stages of each subpopulation of an age-structured SEIRS model with time delay in this paper, which include immature stage and mature stage. In the immature stage, individuals can be born, grow up, die, or survive until the maximum age *a*1, but, in this case, the individuals are not proliferate until the maximum age *a*1. The age *a*<sup>1</sup> is considered to be the maximum age in the immature stage and it is also considered to be the initial age in the mature stage. In the mature stage, individuals have reached maturity (*a* > *a*1) and can grow up, proliferate, die, or survive to the maximum age *A*. Here, we consider susceptible, exposed, infective, and recovered individuals of two age stages in an age-structured SEIRS model. Let *S*(*a*, *t*), *E*(*a*, *t*), *I*(*a*, *t*), and *R*(*a*, *t*) be the distribution densities of susceptible, exposed, infective ,and recovered individuals of age *x* at time *t*, and they take place in the domain Ω = {(*a*, *<sup>t</sup>*) : 0 <sup>≤</sup> *<sup>a</sup>* <sup>≤</sup> *<sup>A</sup>*, 0 <sup>≤</sup> *<sup>t</sup>* <sup>≤</sup> *<sup>T</sup>*}. The integrals *Ns*(*t*) = - *A* <sup>0</sup> *<sup>S</sup>*(*a*, *<sup>t</sup>*)*da*, *Ne*(*t*) = - *A* <sup>0</sup> *E*(*a*, *t*)*da*, *Ni*(*t*) = - *A* <sup>0</sup> *<sup>I</sup>*(*a*, *<sup>t</sup>*)*da*, and *Nr*(*t*) = - *A* <sup>0</sup> *R*(*a*, *t*)*da* are considered as the number of susceptible, exposed, infective and recovered individuals, respectively. The total population *N*(*t*) is: *N*(*t*) = *Ns*(*t*) + *Ne*(*t*) + *Ni*(*t*) + *Nr*(*t*). What needs to be added is that *S*(*a*, *t*), *E*(*a*, *t*), *I*(*a*, *t*), and *R*(*a*, *t*) should belong to *L*1(Ω) because we assume that the initial total population is limited. Then, the following differential equations with time delay are

$$\begin{cases} \frac{\partial S}{\partial t} + \frac{\partial S}{\partial a} = -\left(\pounds(a,t) + \sigma \hat{\partial}\_{1}(a,t)\right) S(a,t) - \left(\gamma\_{1}(a,t-\tau)\int\_{0}^{A} \beta(a,d,t-\tau)I(a',t-\tau) da'\right) S(a,t) \\ \quad + \rho(a,t)R(a,t), (a,t) \in \Omega, \\ \frac{\partial E}{\partial t} + \frac{\partial E}{\partial a} = -\left(\pounds(a,t) + \sigma \hat{\partial}\_{2}(a,t)\right) E(a,t) + \left(\gamma\_{1}(a,t-\tau)\int\_{0}^{A} \beta(a,d,t-\tau)I(a',t-\tau) da'\right) S(a,t) \\ \quad - \gamma\_{2}(a,t)E(a,t), (a,t) \in \Omega, \\ \frac{\partial I}{\partial t} + \frac{\partial I}{\partial a} = -\left(\pounds(a,t) + \sigma \hat{\partial}\_{3}(a,t)\right) I(a,t) + \gamma\_{2}(a,t)E(a,t) - \gamma\_{3}(a,t)I(a,t), (a,t) \in \Omega, \\ \frac{\partial R}{\partial t} + \frac{\partial R}{\partial a} = -\left(\pounds\_{4}(a,t) + \sigma \hat{\partial}\_{4}(a,t)\right) R(a,t) + \gamma\_{3}(a,t)I(a,t) - \rho(a,t)R(a,t), (a,t) \in \Omega, \end{cases} \tag{1}$$

with initial conditions

$$\begin{cases} \mathcal{S}(a,0) = \mathcal{S}\_0(a), a \in [0, A]\_\prime \\ E(a,0) = E\_0(a), a \in [0, A]\_\prime \\ I(a,t) = I\_0(a,t), a \in [0, A]\_\prime, t \in [-\tau, 0]\_\prime \\ R(a,0) = 0, x \in [0, A]\_\prime \end{cases} \tag{2}$$

and boundary conditions

$$\begin{cases} S(0,t) = \int\_{a\_1}^{A} \left( \mu\_1 \theta\_1(a,t) S(a,t) + \mu\_2 \theta\_2(a,t) \left(1 - p(a,t)\right) E(a,t) + \mu\_3 \theta\_3(a,t) \left(1 - q(a,t)\right) I(a,t) \right) \\ \quad + \mu\_4 \theta\_4(a,t) R(a,t) \right) da, t \in (0,T), \\ E(0,t) = \mu\_2 \int\_{a\_1}^{A} p(a,t) \theta\_2(a,t) E(a,t) da, t \in (0,T), \\ I(0,t) = \mu\_3 \int\_{a\_1}^{A} q(a,t) \theta\_3(a,t) I(a,t) da, t \in (0,T), \\ R(0,t) = 0, t \in (0,T), \end{cases} (3)$$

where ˆ *θi*(*a*, *t*)(*i* = 1, 2, 3, 4) denote fertility rates of females of each subpopulation of age *a*; ˆ *θi*(*a*, *t*) = *<sup>θ</sup>i*(*a*, *<sup>t</sup>*), if *<sup>a</sup>* <sup>∈</sup> [*a*1, *<sup>A</sup>*]; <sup>ˆ</sup> *θi*(*a*, *t*) = 0, if *a* ∈/ [*a*1, *A*]; *α*ˆ*i*(*a*, *t*) are natural death rates of each subpopulation of age *a*; *σ* is the provided parameter: *σ* = 1 if the individuals die when they produce and *σ* = 0 if the individuals continue to survive when they produce; *γ*1(*a*, *t*) is a transmission coefficient which describes the varying probability of infectiousness and it is related to a great many social, environmental, and epidemiological factors; *β*(*a*, *a* , *t*) is the contact rate between infected population (age *a* ) and susceptible population (age *a*) per unit time; *τ* (*τ* > 0) is a fixed incubation period of infection; *γ*2(*a*, *t*) is the conversion rate from exposed population to infected population of age *a*; *γ*3(*a*, *t*) is the recovery rate of age *a*; *ρ*(*a*, *t*) is the conversion rate from recovered population losing immunity to the susceptible population of age *a*; *μi*(*i* = 1, 2, 3, 4) are reproductive rates of each subpopulation in the proliferating stage; *p*(*a*, *t*) is part of exposed individuals which procreate new exposed individuals of age *a*; and, similarly, *q*(*a*, *t*) is part of infective individuals which procreate new infective individuals of age *a*.

## **3. Traveling Wave Solution**

In this section, we mainly use the method of characteristics [2,18,25–28] for first-order hyperbolic partial differential equations (Equation (1)). Then, the system in Equation (1) is reduced to nonlinear delayed integro-differential equations along the characteristics curve *a* − *t* = *constant* [25,26]. Let *u* = *a* − *t*, the following system with time delay is obtained:

$$\begin{cases} S\_{l} = -\left(\pounds\_{1}(u+t,t) + \sigma\hat{\theta}\_{1}(u+t,t)\right)S(u+t,t) - \left(\gamma\_{1}(u+t,t-\tau)\int\_{0}^{A}\beta(a,d',t-\tau)I(a',t-\tau)da'\right) \\\\ \times S(u+t,t) + \rho(u+t,t)\mathbb{R}(u+t,t), \\\\ E\_{l} = -\left(\pounds\_{2}(u+t,t) + \sigma\hat{\theta}\_{2}(u+t,t)\right)E(u+t,t) + \left(\gamma\_{1}(u+t,t-\tau)\int\_{0}^{A}\beta(a,d',t-\tau)I(a',t-\tau)da'\right) \\\\ \times S(u+t,t) - \gamma\_{2}(u+t,t)E(u+t,t), \\ I\_{l} = -\left(\pounds\_{3}(u+t,t) + \sigma\hat{\theta}\_{3}(u+t,t)\right)I(u+t,t) + \gamma\_{2}(u+t,t)E(u+t,t) - \gamma\_{3}(u+t,t)I(u+t,t), \\ R\_{l} = -\left(\pounds\_{4}(u+t,t) + \sigma\hat{\theta}\_{4}(u+t,t)\right)R(u+t,t) + \gamma\_{3}(u+t,t)I(u+t,t) - \rho(u+t,t)R(u+t,t), \end{cases} \tag{4}$$

with initial functions

$$\begin{cases} S(\mu, 0) = S\_0(\mu), \\ E(\mu, 0) = E\_0(\mu), \\ I(\mu + t, t) = I\_0(\mu + t, t), t \in [-\tau, 0], \\ R(\mu, 0) = R\_0(\mu) = 0. \end{cases} \tag{5}$$

In general, we divide time interval [0, *<sup>T</sup>*] into *<sup>K</sup>* intervals [*tk*−1, *tk*], where *<sup>k</sup>* = 1, ..., *<sup>K</sup>*, *<sup>t</sup>*<sup>0</sup> = 0, *tK* = *T*, *tk* = *kA*. Then, Ω can be grouped into two sets:

$$
\Omega\_k^{(1)} = \{(a, t) | t \in [(k - 1)A, (a + (k - 1)A)], a \in [0, A]\},\tag{6}
$$

$$
\Omega\_k^{(2)} = \left\{ (a, t) | t \in [(a + (k - 1)A), kA], a \in [0, A] \right\}, \tag{7}
$$

$$
\Omega = \cup\_{k=1}^{K} (\Omega\_k^{(1)} \cup \Omega\_k^{(2)}).\tag{8}
$$

We also need to define an auxiliary set for *k* = 1, ..., *K*:

$$\bar{\Omega}^{(k)} = \left\{ \left[ -\boldsymbol{u}\_n^{(k)}, -\boldsymbol{u}\_{n+1}^{(k)} \right] \mid \boldsymbol{u}\_n^{(k)} = n\boldsymbol{a}\_1 + (k-1)\boldsymbol{A}, \boldsymbol{u} = 1, \ldots, N-1, \boldsymbol{u}\_N^{(k)} = k\boldsymbol{A} \right\},\tag{9}$$

where *N* = [ *<sup>A</sup> a*1 ] + 1, *i f <sup>A</sup> <sup>a</sup>*<sup>1</sup> <sup>−</sup> [ *<sup>A</sup> a*1 ] > 0; *N* = *<sup>A</sup> a*1 , *i f <sup>A</sup> <sup>a</sup>*<sup>1</sup> <sup>−</sup> [ *<sup>A</sup> a*1 ] > 0. [*x*] denotes a whole part of real number *x*. The following hypotheses are given:

(*H*1): *S*0(*a*), *E*0(*a*) and *I*0(*a*, 0) are non-negative and continuous when *a* ∈ [ 0, *A*]; when *a* → *A* − 0, *S*0(*A*) = 0, *E*0(*A*) = 0 and *I*0(*A*, 0) = 0.

(*H*2): *α*ˆ*i*(*a*, *t*), *θi*(*a*, *t*), *γi*(*a*, *t*), *β*(*a*, *t*), *ρ*(*a*, *t*), *p*(*a*, *t*), *q*(*a*, *t*) ∈ *C*(Ω).

(*H*3): *<sup>α</sup>*ˆ*i*(*a*, *<sup>t</sup>*) <sup>&</sup>gt; 0, *<sup>θ</sup>i*(*a*, *<sup>t</sup>*) <sup>≥</sup> 0, *<sup>γ</sup>i*(*a*, *<sup>t</sup>*) <sup>&</sup>gt; 0, *<sup>β</sup>*(*a*, *<sup>t</sup>*) <sup>≥</sup> 0. In addition, 0 <sup>&</sup>lt; - *A <sup>a</sup>*<sup>1</sup> *<sup>θ</sup>i*(*a*, *<sup>t</sup>*)*da* ≤ 1, 0 < - *A* <sup>0</sup> *β*(*a*, *t*)*dx* ≤ 1, 0 ≤ *ρ*(*a*, *t*) ≤ 1, 0 ≤ *p*(*a*, *t*) ≤ 1, 0 ≤ *q*(*a*, *t*) ≤ 1. (*H*4): zero-order compatibility conditions are

$$\begin{cases} S\_0(0) = \int\_{a\_1}^A \left( \mu\_1 \theta\_1(a,0)S\_0(a) + \mu\_2 \theta\_2(a,0) \left(1 - p(a,0)\right) E\_0(a) + \mu\_3 \theta\_3(a,0) \left(1 - q(a,0)\right) I\_0(a,0) \right) da, \\ E\_0(0) = \mu\_2 \int\_{a\_1}^A \theta\_2(a,0) p(a,0) E\_0(a) da, \\ I\_0(0,0) = \mu\_3 \int\_{a\_1}^A \theta\_3(a,0) q(a,0) I\_0(a,0) da, \\ R\_0(0) = 0. \end{cases} \tag{10}$$

For convenience, we assume Hypotheses (*H*1)–(*H*4) are satisfied.

Then, we solve the system in Equation (4) by using the well-known steps method [8,9]. For the first step *h* = 1, i.e., for the first time interval *t* ∈ [0, *τ*], the solution to the initial value problem in Equation (4) can be obtained according to the known values of function *I*(*x*, *t* − *τ*). Repeating the pervious step for *h* = 2, 3, 4, ..., *H*, for the time interval *t* ∈ [(*h* − 1)*τ*, *hτ*], we can obtain the solution to the initial value problem in Equation (4) for the whole time interval *t* ∈ [0, *T*]:

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *<sup>S</sup>*1(*u*, *<sup>t</sup>*) = *<sup>S</sup>*0(*u*)*W*1(*u*, 0, *<sup>t</sup>*, *<sup>τ</sup>*) + *<sup>t</sup>* 0 *W*1(*u*, *φ*, *t*, *τ*)*ρ*(*u* + *φ*, *φ*)*R*(*u* + *φ*, *φ*)*dφ*, *t* ∈ [0, *τ*], *Sh*(*u*, *<sup>t</sup>*) = *Sh*−1(*u*)*W*1(*u*,(*<sup>h</sup>* <sup>−</sup> <sup>1</sup>)*τ*, *<sup>t</sup>*, *<sup>τ</sup>*) + *<sup>t</sup>* (*h*−1)*τ W*1(*u*, *φ*, *t*, *τ*)*ρ*(*u* + *φ*, *φ*)*R*(*u* + *φ*, *φ*)*dφ* <sup>=</sup> *<sup>S</sup>*0(*u*)*W*1(*v*, 0, *<sup>t</sup>*, *<sup>τ</sup>*) + *<sup>t</sup>* 0 *W*1(*u*, *φ*, *t*, *τ*)*ρ*(*v* + *φ*, *φ*)*R*(*u* + *φ*, *φ*)*dφ*, *t* ∈ [(*h* − 1)*τ*, *hτ*], *<sup>W</sup>*1(*u*, *<sup>t</sup>*0, *<sup>t</sup>*, *<sup>τ</sup>*) = *exp* − *t t*0 *α*1(*u* + *φ*, *φ*) + *γ*1(*u* + *φ*, *φ* − *τ*)*D*(*u* + *φ*, *φ* − *τ*) *dφ* , (11) ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *<sup>E</sup>*1(*u*, *<sup>t</sup>*) = *<sup>E</sup>*0(*u*)*W*2(*u*, 0, *<sup>t</sup>*) + *<sup>t</sup>* 0 *W*2(*u*, *φ*, *t*)*γ*1(*u* + *φ*, *φ* − *τ*)*S*(*u* + *φ*, *φ*)*D*(*u* + *φ*, *φ* − *τ*)*dφ*, *t* ∈ [0, *τ*], *Eh*(*u*, *<sup>t</sup>*) = *Eh*−1(*u*)*W*2(*u*,(*<sup>h</sup>* <sup>−</sup> <sup>1</sup>)*τ*, *<sup>t</sup>*) + *<sup>t</sup>* (*h*−1)*τ W*2(*u*, *φ*, *t*)*γ*1(*u* + *φ*, *φ* − *τ*)*S*(*u* + *φ*, *φ*) <sup>×</sup> *<sup>D</sup>*(*<sup>u</sup>* <sup>+</sup> *<sup>φ</sup>*, *<sup>φ</sup>* <sup>−</sup> *<sup>τ</sup>*)*d<sup>φ</sup>* <sup>=</sup> *<sup>E</sup>*0(*u*)*W*2(*u*, 0, *<sup>t</sup>*) + *<sup>t</sup>* 0 *W*2(*u*, *φ*, *t*)*γ*1(*u* + *φ*, *φ* − *τ*) × *S*(*u* + *φ*, *φ*)*D*(*u* + *φ*, *φ* − *τ*)*dφ*, *t* ∈ [(*h* − 1)*τ*, *hτ*], *<sup>W</sup>*2(*u*, *<sup>t</sup>*0, *<sup>t</sup>*) = *exp* − *t t*0 *α*2(*u* + *φ*, *φ*)+*γ*2(*u* + *φ*, *φ*) *dφ* , (12) ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *<sup>I</sup>*1(*u*, *<sup>t</sup>*) = *<sup>I</sup>*0(*u*, 0)*W*3(*u*, 0, *<sup>t</sup>*) + *<sup>t</sup>* 0 *W*3(*u*, *φ*, *t*)*γ*2(*u* + *φ*, *φ* − *τ*)*E*(*u* + *φ*, *φ*)*dφ*, *t* ∈ [0, *τ*], *Ih*(*u*, *<sup>t</sup>*) = *Ih*−1(*u*,(*<sup>h</sup>* <sup>−</sup> <sup>1</sup>)*τ*)*W*3(*u*,(*<sup>h</sup>* <sup>−</sup> <sup>1</sup>)*τ*, *<sup>t</sup>*) + *<sup>t</sup>* (*h*−1)*τ W*3(*u*, *φ*, *t*)*γ*2(*u* + *φ*, *φ* − *τ*) <sup>×</sup> *<sup>E</sup>*(*<sup>u</sup>* <sup>+</sup> *<sup>φ</sup>*, *<sup>φ</sup>*)*d<sup>φ</sup>* <sup>=</sup> *<sup>I</sup>*0(*u*, 0)*W*3(*u*, 0, *<sup>t</sup>*) + *<sup>t</sup>* 0 *W*3(*u*, *φ*, *t*)*γ*2(*u* + *φ*, *φ* − *τ*) × *E*(*u* + *φ*, *φ*)*dφ*, *t* ∈ [(*h* − 1)*τ*, *hτ*], *<sup>W</sup>*3(*u*, *<sup>t</sup>*0, *<sup>t</sup>*) = *exp* − *t t*0 *α*3(*u* + *φ*, *φ*)+*γ*3(*u* + *φ*, *φ*) *dφ* , (13) ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *<sup>R</sup>*1(*u*, *<sup>t</sup>*) = *<sup>t</sup>* 0 *W*4(*u*, *φ*, *t*)*γ*3(*u* + *φ*, *φ* − *τ*)*I*(*u* + *φ*, *φ*)*dφ*, *t* ∈ [0, *τ*] *Rh*(*u*, *<sup>t</sup>*) = *Rh*−1(*u*,(*<sup>h</sup>* <sup>−</sup> <sup>1</sup>)*τ*)*W*4(*u*,(*<sup>h</sup>* <sup>−</sup> <sup>1</sup>)*τ*), *<sup>t</sup>*) + *<sup>t</sup>* (*h*−1)*τ W*4(*u*, *φ*, *t*)*γ*3(*u* + *φ*, *φ* − *τ*) × *I*(*u* + *φ*, *φ*)*dφ* = *t* 0 *W*4(*u*, *φ*, *t*)*γ*3(*u* + *φ*, *φ* − *τ*)*I*(*u* + *φ*, *φ*)*dφ*, *t* ∈ [(*h* − 1)*τ*, *hτ*], *<sup>W</sup>*4(*u*, *<sup>t</sup>*0, *<sup>t</sup>*) = *exp* − *t t*0 *α*4(*u* + *φ*, *φ*) + *ρ*(*u* + *φ*, *φ*) *dφ* , (14)

where *αi*(*a*, *t*) = *α*ˆ*i*(*a*, *t*) + *σ* ˆ *<sup>θ</sup>i*(*a*, *<sup>t</sup>*), *<sup>i</sup>* <sup>=</sup> 1, 2, 3, 4, *<sup>D</sup>*(*<sup>u</sup>* <sup>+</sup> *<sup>φ</sup>*, *<sup>φ</sup>* <sup>−</sup> *<sup>τ</sup>*) = - *A* <sup>0</sup> *β*(*u* + *φ*, *η*, *φ* − *τ*)*I*(*η*, *φ* − *τ*)*dη*.

The exact solution to Equations (11)–(14) of the system in Equation (4) can be expressed in terms of the original variables. Let *k* = 1 and *u* = *a* − *t*; we have the solution of the system in Equation (1) in the sets Ω(1) *<sup>k</sup>* and <sup>Ω</sup>(2) *<sup>k</sup>* . Repeating this process for *k* = 2, 3, ...*K*, the form of explicit traveling wave solution, i.e., numerical solution of the system in Equation (1), can be obtained in all domains Ω(1) *<sup>k</sup>* and Ω(2) *<sup>k</sup>* for the values of *k* = 1, ..., *K*:

$$(S(a,t),E(a,t),I(a,t),R(a,t)) = \begin{cases} S\_h^{(k-1)}(u,t) = S\_h^{(k-1)}(a-t,t), \\ E\_h^{(k-1)}(u,t) = E\_h^{(k-1)}(a-t,t), \\ I\_h^{(k-1)}(u,t) = I\_h^{(k-1)}(a-t,t), \\ R\_h^{(k-1)}(u,t) = R\_h^{(k-1)}(a-t,t), \\ S\_h^{(k)}(u,t) = S\_h^{(k)}(a-t,t), \\ E\_h^{(k)}(u,t) = E\_h^{(k)}(a-t,t), \\ I\_h^{(k)}(u,t) = I\_h^{(k)}(a-t,t), \\ R\_h^{(k)}(u,t) = R\_h^{(k)}(a-t,t), \\ \end{cases} \tag{15}$$

$$\begin{cases} \mathcal{S}^{(0)}(\boldsymbol{u},t) = \mathcal{S}\_{0}(\boldsymbol{u})\mathcal{W}\_{1}(\boldsymbol{u},0,t,\tau) + \int\_{0}^{t} \mathcal{W}\_{1}(\boldsymbol{u},\boldsymbol{\phi},t,\tau)\rho(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})\mathcal{R}^{(0)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})d\boldsymbol{\phi}, \\\mathcal{S}^{(k)}(\boldsymbol{u},t) = \mathcal{F}^{(k)}(\boldsymbol{u})\mathcal{W}\_{1}(\boldsymbol{u},-\boldsymbol{u},t,\tau) + \int\_{-\boldsymbol{u}}^{t} \mathcal{W}\_{1}(\boldsymbol{u},\boldsymbol{\phi},t,\tau)\rho(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})\mathcal{R}^{(k)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})d\boldsymbol{\phi}, \end{cases} \tag{16}$$

$$\begin{cases} E^{(0)}(\boldsymbol{u},t) = \mathbb{E}\_{0}(\boldsymbol{u})\mathcal{W}\_{2}(\boldsymbol{u},0,t) + \int\_{0}^{t} \mathcal{W}\_{2}(\boldsymbol{u},\boldsymbol{\phi},t)\gamma\_{1}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau})\mathcal{S}^{(0)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})\mathcal{D}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau})d\boldsymbol{\phi}, \\\ E^{(k)}(\boldsymbol{u},t) = \mathcal{G}^{(k)}(\boldsymbol{u})\mathcal{W}\_{2}(\boldsymbol{u},-\boldsymbol{u},t) + \int\_{-\boldsymbol{u}}^{t} \mathcal{W}\_{2}(\boldsymbol{u},\boldsymbol{\phi},t)\gamma\_{1}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau})\mathcal{S}^{(k)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})\mathcal{D}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau})d\boldsymbol{\phi}, \end{cases} \tag{17}$$

$$\begin{cases} I^{(0)}(\boldsymbol{u},t) = l\_0(\boldsymbol{u},0)\mathcal{W}\_3(\boldsymbol{u},0,t) + \int\_0^t \mathcal{W}\_3(\boldsymbol{u},\boldsymbol{\phi},t)\gamma\_2(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau})E^{(0)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})d\boldsymbol{\phi}, \\\\ I^{(k)}(\boldsymbol{u},t) = Q^{(k)}(\boldsymbol{u})\mathcal{W}\_3(\boldsymbol{u},-\boldsymbol{u},t) + \int\_{-\boldsymbol{u}}^t \mathcal{W}\_3(\boldsymbol{u},\boldsymbol{\phi},t)\gamma\_2(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau})E^{(k)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})d\boldsymbol{\phi}, \\\\ I^{(0)}(\boldsymbol{u},\boldsymbol{\omega}) = \int\_{-\boldsymbol{u}}^t \mathcal{W}\_3(\boldsymbol{u},\boldsymbol{\omega},\boldsymbol{\phi},t)\gamma\_2(\boldsymbol{u},\boldsymbol{\phi},\boldsymbol{\phi},\boldsymbol{\phi})\gamma\_2(\boldsymbol{u},\boldsymbol{\phi},\boldsymbol{\phi},\boldsymbol{\phi})d\boldsymbol{\phi}, \end{cases} \tag{18}$$

$$\begin{cases} \mathcal{R}^{(0)}(\boldsymbol{u},t) = \int\_{0}^{t} \mathcal{W}\_{4}(\boldsymbol{u},\boldsymbol{\phi},t) \gamma\_{3}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau}) I^{(0)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}) d\boldsymbol{\phi}, \\\mathcal{R}^{(k)}(\boldsymbol{u},t) = \int\_{-\boldsymbol{u}}^{t} \mathcal{W}\_{4}(\boldsymbol{u},\boldsymbol{\phi},t) \gamma\_{3}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau}) I^{(k)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}) d\boldsymbol{\phi}, \end{cases} \tag{19}$$

where *W*1(*u*, *t*0, *t*, *τ*), *W*2(*u*, *t*0, *t*), *W*3(*u*, *t*0, *t*), and *W*4(*u*, *t*0, *t*) are shown in Equations (11)–(14), respectively. *<sup>F</sup>*(*k*)(*u*), *<sup>G</sup>*(*k*)(*u*), and *<sup>Q</sup>*(*k*)(*u*) are given by defining functions *<sup>F</sup>*(*k*) *<sup>n</sup>* (*u*), *<sup>G</sup>*(*k*) *<sup>n</sup>* (*u*), *<sup>Q</sup>*(*k*) *<sup>n</sup>* (*u*), *S*(*k*) *<sup>n</sup>* (*u*), *<sup>E</sup>*(*k*) *<sup>n</sup>* (*u*), *I* (*k*) *<sup>n</sup>* (*u*), and *<sup>R</sup>*(*k*) *<sup>n</sup>* (*u*), *k* = 1, ..., *K*:

$$(F^{(k)}(\mathfrak{u})), G^{(k)}(\mathfrak{u}), Q^{(k)}(\mathfrak{u})) = (F\_{\mathfrak{n}}^{(k)}(\mathfrak{u})), G\_{\mathfrak{n}}^{(k)}(\mathfrak{u}), Q\_{\mathfrak{n}}^{(k)}(\mathfrak{u})), \mathfrak{u} \in [-\mathfrak{u}\_{\mathfrak{n}}^{(k)}, -\mathfrak{u}\_{\mathfrak{n}-1}^{(k)}], \mathfrak{u} = 1, \dots, N\_{\mathsf{T}}$$

thus we have

$$\begin{cases} \mathcal{S}\_{n}^{(k)}(\boldsymbol{u},t) = F\_{n}^{(k)}(\boldsymbol{u})\mathcal{W}\_{1}(\boldsymbol{u},-\boldsymbol{u},t,\tau) + \int\_{-\mathfrak{u}}^{t} \mathcal{W}\_{1}(\boldsymbol{u},\boldsymbol{\phi},t,\tau)\boldsymbol{\rho}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})\mathcal{R}\_{n}^{(k)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})d\boldsymbol{\phi}, \\\ E\_{n}^{(k)}(\boldsymbol{u},t) = G\_{n}^{(k)}(\boldsymbol{u})\mathcal{W}\_{2}(\boldsymbol{u},-\boldsymbol{u},t) + \int\_{-\mathfrak{u}}^{t} \mathcal{W}\_{2}(\boldsymbol{u},\boldsymbol{\phi},t)\gamma\_{1}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau})\mathcal{S}\_{n}^{(k)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})D(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau})d\boldsymbol{\phi}, \\\ I\_{n}^{(k)}(\boldsymbol{u},t) = Q\_{n}^{(k)}(\boldsymbol{u})\mathcal{W}\_{3}(\boldsymbol{u},-\boldsymbol{u},t) + \int\_{-\mathfrak{u}}^{t} \mathcal{W}\_{3}(\boldsymbol{u},\boldsymbol{\phi},t)\gamma\_{2}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau})\mathcal{E}\_{n}^{(k)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})d\boldsymbol{\phi}, \\\ R\_{n}^{(k)}(\boldsymbol{u},t) = \int\_{-\mathfrak{u}}^{t} \mathcal{W}\_{4}(\boldsymbol{u},\boldsymbol{\phi},t)\gamma\_{3}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi}-\boldsymbol{\tau})\mathcal{I}\_{n}^{(k)}(\boldsymbol{u}+\boldsymbol{\phi},\boldsymbol{\phi})d\boldsymbol{\phi}, \end{cases}$$

where *<sup>n</sup>* <sup>=</sup> 1, ..., *<sup>N</sup><sup>τ</sup>* <sup>−</sup> 1 and functions *<sup>F</sup>*(*k*) *<sup>n</sup>* (*u*), *<sup>G</sup>*(*k*) *<sup>n</sup>* (*u*), and *<sup>Q</sup>*(*k*) *<sup>n</sup>* (*u*) can be defined according to the recurrent algorithm as follows:

*F*(*k*) <sup>1</sup> (*v*) = *<sup>A</sup>*+*<sup>v</sup> a*1+*v <sup>μ</sup>*1*θ*1(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*S*(*k*−1) (*u*, <sup>−</sup>*v*) + *<sup>μ</sup>*2*θ*2(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)(<sup>1</sup> <sup>−</sup> *<sup>p</sup>*(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)) *<sup>E</sup>*(*k*−1) (*u*, −*v*) + *μ*3*θ*3(*u* − *v*, −*v*)(1 − *q*(*u* − *v*, −*v*)) *I* (*k*−1) (*u*, <sup>−</sup>*v*) + *<sup>μ</sup>*4*θ*4(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*R*(*k*−1) (*u*, −*v*) *du*, *<sup>v</sup>* <sup>∈</sup> [−*u*(*k*) <sup>1</sup> , <sup>−</sup>*u*(*k*) <sup>0</sup> ], *G*(*k*) <sup>1</sup> (*v*) = *μ*<sup>2</sup> *A*+*v a*1+*v <sup>p</sup>*(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*θ*2(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*u*)*E*(*k*−1) (*u*, <sup>−</sup>*v*)*du*, *<sup>v</sup>* <sup>∈</sup> [−*u*(*k*) <sup>1</sup> , <sup>−</sup>*u*(*k*) <sup>0</sup> ], *Q*(*k*) <sup>1</sup> (*v*) = *μ*<sup>3</sup> *A*+*v a*1+*v q*(*u* − *v*, −*v*)*θ*3(*u* − *v*, −*u*)*I* (*k*−1) (*u*, <sup>−</sup>*v*)*du*, *<sup>v</sup>* <sup>∈</sup> [−*u*(*k*) <sup>1</sup> , <sup>−</sup>*u*(*k*) <sup>0</sup> ], *F*(*k*) <sup>2</sup> (*v*) = <sup>−</sup>*u*(*k*) 0 *a*1+*v <sup>μ</sup>*1*θ*1(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*S*(*k*) (*u*, <sup>−</sup>*v*) + *<sup>μ</sup>*2*θ*2(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)(<sup>1</sup> <sup>−</sup> *<sup>p</sup>*(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)) *<sup>E</sup>*(*k*) (*u*, −*v*) + *μ*3*θ*3(*u* − *v*, −*v*)(1 − *q*(*u* − *v*, −*v*)) *I* (*k*) (*u*, <sup>−</sup>*v*) + *<sup>μ</sup>*4*θ*4(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*R*(*k*) (*u*, −*v*) *du* + *A*+*v* −*u*(*k*) 0 *<sup>μ</sup>*1*θ*1(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*S*(*k*−1) (*u*, −*v*) + *μ*2*θ*2(*u* − *v*, −*v*)(1 − *p*(*u* − *v*, −*v*)) <sup>×</sup> *<sup>E</sup>*(*k*−1) (*u*, −*v*) + *μ*3*θ*3(*u* − *v*, −*v*)(1 − *q*(*u* − *v*, −*v*)) *I* (*k*−1) (*u*, −*v*) + *μ*4*θ*4(*u* − *v*, −*v*) <sup>×</sup> *<sup>R</sup>*(*k*−1) (*u*, −*v*) *du*, *<sup>v</sup>* <sup>∈</sup> [−*u*(*k*) <sup>2</sup> , <sup>−</sup>*u*(*k*) <sup>1</sup> ], *G*(*k*) <sup>2</sup> (*v*) = *μ*<sup>2</sup> <sup>−</sup>*u*(*k*) 0 *a*1+*v <sup>p</sup>*(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*θ*2(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*u*)*E*(*k*) (*u*, −*v*)*du* + *μ*<sup>2</sup> *A*+*v* −*u*(*k*) 0 *p*(*u* − *v*, −*v*) <sup>×</sup> *<sup>θ</sup>*2(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*u*)*E*(*k*−1) (*u*, <sup>−</sup>*v*)*du*, *<sup>v</sup>* <sup>∈</sup> [−*u*(*k*) <sup>2</sup> , <sup>−</sup>*u*(*k*) <sup>1</sup> ], *Q*(*k*) <sup>2</sup> (*v*) = *μ*<sup>3</sup> <sup>−</sup>*u*(*k*) 0 *a*1+*v q*(*u* − *v*, −*v*)*θ*3(*u* − *v*, −*u*)*I* (*k*) (*u*, −*v*)*du* + *μ*<sup>3</sup> *A*+*v* −*u*(*k*) 0 *q*(*u* − *v*, −*v*) × *θ*3(*u* − *v*, −*u*)*I* (*k*−1) (*u*, <sup>−</sup>*v*)*du*, *<sup>v</sup>* <sup>∈</sup> [−*u*(*k*) <sup>2</sup> , <sup>−</sup>*u*(*k*) <sup>1</sup> ], *<sup>F</sup>*(*k*) *<sup>n</sup>* (*v*) = <sup>−</sup>*u*(*k*) *n*−2 *a*1+*v <sup>μ</sup>*1*θ*1(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*S*(*k*) (*u*, <sup>−</sup>*v*) + *<sup>μ</sup>*2*θ*2(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)(<sup>1</sup> <sup>−</sup> *<sup>p</sup>*(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)) *<sup>E</sup>*(*k*) (*u*, −*v*) + *μ*3*θ*3(*u* − *v*, −*v*)(1 − *q*(*u* − *v*, −*v*)) *I* (*k*) (*u*, <sup>−</sup>*v*) + *<sup>μ</sup>*4*θ*4(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*R*(*k*) (*u*, −*v*) *du* + *n*−3 ∑ *i*=0 <sup>−</sup>*u*(*k*) *i* −*u*(*k*) *i*+1 *<sup>μ</sup>*1*θ*1(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*S*(*k*) *<sup>i</sup>*+1(*u*, −*v*) + *μ*2*θ*2(*u* − *v*, −*v*)(1 − *p*(*u* − *v*, −*v*)) <sup>×</sup> *<sup>E</sup>*(*k*) *<sup>i</sup>*+1(*u*, −*v*) + *μ*3*θ*3(*u* − *v*, −*v*)(1 − *q*(*u* − *v*, −*v*)) *I* (*k*) *<sup>i</sup>*+1(*u*, −*v*) + *μ*4*θ*4(*u* − *v*, −*v*) + *R*(*k*) *<sup>i</sup>*+1(*u*, −*v*) *du* + *A*+*v* −*u*(*k*) 0 *<sup>μ</sup>*1*θ*1(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*S*(*k*−1) (*u*, −*v*) + *μ*2*θ*2(*u* − *v*, −*v*) <sup>×</sup> (<sup>1</sup> <sup>−</sup> *<sup>p</sup>*(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)) *<sup>E</sup>*(*k*−1) (*u*, −*v*) + *μ*3*θ*3(*u* − *v*, −*v*)(1 − *q*(*u* − *v*, −*v*)) *I* (*k*−1) (*u*, −*v*) <sup>+</sup> *<sup>μ</sup>*4*θ*4(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*R*(*k*−1) (*u*, −*v*) *du*, *<sup>v</sup>* <sup>∈</sup> [−*u*(*k*) *<sup>n</sup>* , <sup>−</sup>*u*(*k*) *<sup>n</sup>*−1], *<sup>n</sup>* <sup>=</sup> 3, ..., *<sup>N</sup>τ*, *G*(*k*) *<sup>n</sup>* (*v*) = *μ*<sup>2</sup> <sup>−</sup>*u*(*k*) *n*−2 *a*1+*v <sup>p</sup>*(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*θ*2(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*u*)*E*(*k*) *<sup>n</sup>*−1(*u*, <sup>−</sup>*v*)*du* <sup>+</sup> *<sup>μ</sup>*<sup>2</sup> *n*−3 ∑ *i*=0 <sup>−</sup>*u*(*k*) *i* −*u*(*k*) *i*+1 *p*(*u* − *v*, −*v*) <sup>×</sup> *<sup>θ</sup>*2(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*u*)*E*(*k*) *<sup>i</sup>*+1(*u*, −*v*)*du* + *μ*<sup>2</sup> *A*+*v* −*u*(*k*) 0 *<sup>p</sup>*(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*v*)*θ*2(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*, <sup>−</sup>*u*)*E*(*k*−1) (*u*, −*v*)*du*, *<sup>v</sup>* <sup>∈</sup> [−*u*(*k*) *<sup>n</sup>* , <sup>−</sup>*u*(*k*) *<sup>n</sup>*−1], *<sup>n</sup>* <sup>=</sup> 3, ..., *<sup>N</sup>τ*,

*Mathematics* **2020**, *8*, 455

$$\begin{split} &Q\_{n}^{(k)}(v) = \mu\_{3} \int\_{u\_{1}+v}^{-u\_{2}^{(k)}} q(u-v,-v) \theta\_{3}(u-v,-u) I\_{n-1}^{(k)}(u\_{\prime}-v) du + \mu\_{3} \sum\_{i=0}^{n-3} \int\_{-u\_{i+1}^{(k)}}^{-u\_{i}^{(k)}} q(u-v,-v) \\ &\times \theta\_{3}(u-v,-u) I\_{i+1}^{(k)}(u\_{\prime}-v) du + \mu\_{3} \int\_{-u\_{0}^{(k)}}^{A+v} q(u-v,-v) \theta\_{3}(u-v,-u) I^{(k-1)}(u\_{\prime}-v) du, \\ &v \in [-u\_{n}^{(k)},-u\_{n-1}^{(k)}], n = 3, \dots, N\_{\tau}, \\ &F\_{1}^{(k)}(u) = \int\_{u\_{1}}^{A} \left[ \mu\_{1} \theta\_{1}(v,-u) S^{(k-1)}(v,-u) + \mu\_{2} \theta\_{2}(v,-u) \left(1 - p(v,-u)\right) E^{(k-1)}(v,-u) \right. \\ &+ \mu\_{3} \theta\_{3}(v,-u) \left(1 - q(v,-u)\right) I^{(k-1)}(v,-u) + \mu\_{4} \theta\_{4}(v,-u) R^{(k-1)}(v,-u) \Big] dv, \\ &G^{(k)}(u) = \mu\_{2} \int\_{u\_{1}}^{A} p(v,-u) \theta\_{2}(v,-u) E^{(k-1)}(v,-u) dv, \\ &Q^{(k)}(u) = \mu\_{3} \int\_{u\_{1}}^{A} q(v,-u) \theta\_{3}(v,-u) I^{(k-1)}(v,-u) dv. \end{split}$$

**Theorem 1.** *If Assumptions* (*H*1)*–*(*H*4) *are satisfied, there exists a unique continuous traveling wave solution to Equations (15)–(19) of the system in Equation (1).*

**Proof.** By analogy with [28], if *S*0(*a*), *E*0(*a*), *I*0(*a*, 0), and *R*0(*a*) satisfy compatibility conditions (i.e., Hypothesis *H*4), then two parts of the traveling wave in Equation (15) can be combined continuously. If the parameters of the system in Equation (1) fulfill Hypotheses (*H*1)–(*H*3), there exists a unique continuous traveling wave solution to Equations (15)–(19) of the system in Equation (1). The method is very similar, thus we omit it.

#### **4. Stability Analysis of System**

Consider the nonlinear autonomous system in Equation (1) where the following parameters are constants: *θi*(*a*, *t*) = *θi*0(*i* = 1, 2, 3, 4), *αi*(*a*, *t*) = *αi*<sup>0</sup> = *α<sup>i</sup>* + *σθi*0(*i* = 1, 2, 3, 4), *γ*1(*a*, *t*) = *γ*10, *γ*2(*a*, *t*) = *γ*20, *γ*3(*a*, *t*) = *γ*30, *β*(*a*, *a* , *t*) = *β*0, *p*(*a*, *t*) = *p*0, *q*(*a*, *t*) = *q*0, and *ρ*(*a*, *t*) = *ρ*0, where *αi*0, *θi*0, *γi*0, *β*0, *p*0, *q*0, and *ρ*<sup>0</sup> are positive constants. In this paper, a partial differential equation (Equation (1)) with the initial-boundary values in Equations (2) and (3) reduced to a nonlinear ordinary differential equation. Take the maturity age *a*<sup>1</sup> → 0, which is not an essential simplification. Integrating Equation (1) in regard to age *a* from 0 to *A*, and using the real conditions *S*(*A*, *t*) = 0, *E*(*A*, *t*) = 0, *I*(*A*, *t*) = 0 and *R*(*A*, *t*) = 0, the initial value problem of the nonlinear ordinary differential equation autonomous system that describes the population dynamics of the number of population *Ns*(*t*), *Ne*(*t*), *Ni*(*t*), and *Nr*(*t*) is:

$$\begin{cases} N\_s'(t) = -(\mathfrak{a}\_{10} - \mu\_1 \theta\_{10}) N\_s(t) - \gamma\_{10} \theta\_0 N\_s(t) N\_i(t - \tau) + \mu\_2 \theta\_{20} (1 - p\_0) N\_\varepsilon(t) \\ \quad + \mu\_3 \theta\_{30} (1 - q\_0) N\_i(t) + (\mu\_4 \theta\_{40} + \rho\_0) N\_r(t), \\ N\_\varepsilon'(t) = -(\mathfrak{a}\_{20} + \gamma\_{20} - \mu\_2 \theta\_{20} p\_0) N\_\varepsilon(t) + \gamma\_{10} \theta\_0 N\_\varepsilon(t) N\_i(t - \tau), \\ N\_i'(t) = -(\mathfrak{a}\_{30} + \gamma\_{30} - \mu\_3 \theta\_{30} q\_0) N\_i(t) + \gamma\_{20} N\_\varepsilon(t), \\ N\_r'(t) = -(\mathfrak{a}\_{40} + \rho\_0) N\_r(t) + \gamma\_{30} N\_i(t), \end{cases} \tag{20}$$

with initial functions

$$\begin{cases} N\_s(0) = \int\_0^A S\_0(a) da, \\ N\_\epsilon(t) = \int\_0^A E\_0(a, t) da, t \in [-\tau, 0], \\ N\_i(0) = \int\_0^A I\_0(a) da, \\ N\_r(0) = 0. \end{cases} \tag{21}$$

The basic reproduction number [31] is defined as

$$R\_0 = \frac{N\gamma\_{20}\gamma\_{10}\beta\_0}{(a\_{20} + \gamma\_{20} - \mu\_2\theta\_{20}p\_0)(a\_{30} + \gamma\_{30} - \mu\_3\theta\_{30}q\_0)}.\tag{22}$$

According to the analysis of [28], if *R*<sup>0</sup> > 1, it can lead to the outbreak of infectious diseases. Hence, the following theorem is obtained.

**Theorem 2.** *If γ*<sup>10</sup> = 0*, the system in Equation (20) has only a disease-free equilibrium point H*0*=*(0*,* 0*,* 0*,* 0)*; if R*<sup>0</sup> > 1 *and the parameters of the system in Equation (20) satisfy R*<sup>1</sup> > 1*, R*<sup>2</sup> < 1*, R*<sup>3</sup> < 1*, and R*<sup>4</sup> < 1*, where <sup>R</sup>*<sup>1</sup> = *<sup>μ</sup>*1*θ*<sup>10</sup> *<sup>α</sup>*<sup>10</sup> *, <sup>R</sup>*<sup>2</sup> <sup>=</sup> *<sup>μ</sup>*2*θ*<sup>20</sup> *<sup>p</sup>*<sup>0</sup> *<sup>α</sup>*20+*γ*<sup>20</sup> *, <sup>R</sup>*<sup>3</sup> <sup>=</sup> *<sup>μ</sup>*3*θ*30*q*<sup>0</sup> *<sup>α</sup>*30+*γ*<sup>30</sup> *, and <sup>R</sup>*<sup>4</sup> <sup>=</sup> *<sup>μ</sup>*3*θ*30(1−*q*0)*γ*20+*γ*20*γ*30(*μ*4*θ*40+*ρ*0)(*α*40+*ρ*0)−<sup>1</sup> (*α*20+*γ*20−*μ*2*θ*20)(*α*30+*γ*30−*μ*3*θ*30*q*0) *. The system in Equation (20) has only a positive endemic equilibrium point H*∗ = (*N*∗ *<sup>s</sup> , N*<sup>∗</sup> *<sup>e</sup> , N*<sup>∗</sup> *<sup>i</sup> , N*<sup>∗</sup> *<sup>r</sup>* )*, where N*<sup>∗</sup> *<sup>s</sup> , N*<sup>∗</sup> *<sup>e</sup> , N*<sup>∗</sup> *<sup>i</sup> , and N*∗ *<sup>r</sup> are given in the proof.*

**Proof.** When *γ*<sup>10</sup> = 0, that is no infectious diseases, there exists one disease-free equilibrium point *H*0=(0, 0, 0, 0). When the fertility rate and the death rate of susceptible population are satisfied by

$$R\_1 = \frac{\mu\_1 \theta\_{10}}{\alpha\_{10}} > 1,\tag{23}$$

where *R*<sup>1</sup> is a dimensionless index for the existence of the disease-free equilibrium point. It can be seen that *R*<sup>1</sup> > 1 presents that death rate of susceptible population is less than their reproductive rate. We can take it further into consideration with the dynamic behavior of susceptible population. The endemic equilibrium point *H*∗= (*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) is the solution of nonlinear system

$$- \left(\mu\_{10} - \mu\_1 \theta\_{10}\right) \mathcal{N}\_t^\* - \gamma\_{10} \theta\_0 \mathcal{N}\_s^\* \mathcal{N}\_l^\* + \mu\_2 \theta\_{20} (1 - p\_0) \mathcal{N}\_t^\* + \mu\_3 \theta\_{30} (1 - q\_0) \mathcal{N}\_l(t) + (\mu\_4 \theta\_{40} + \rho\_0) \mathcal{N}\_t^\* = 0,\tag{24}$$

$$- \left( \mu \varepsilon \mathfrak{o} + \gamma \mathfrak{z} \mathfrak{o} - \mu \varepsilon \mathfrak{e} \mathfrak{z} \wp p\_0 \right) \mathcal{N}\_\varepsilon^\* + \gamma\_{10} \mathfrak{z} \wp \mathcal{N}\_\mathfrak{s}^\* \mathcal{N}\_\mathfrak{l}^\* = 0,\tag{25}$$

$$- \left( \mu\_{30} + \gamma\_{30} - \mu\_3 \theta\_{30} \rho\_0 \right) N\_l^\* + \gamma\_{20} N\_\varepsilon^\* = 0,\tag{26}$$

$$- \left( a\_{40} + \rho\_0 \right) N\_r^\* + \gamma\_{50} N\_l^\* = 0. \tag{27}$$

On the basis of the parameters of the system in Equation (1) satisfying Hypothesis (*H*1)–(*H*3), the endemic equilibrium point *N*∗ *<sup>s</sup>* > 0, *N*<sup>∗</sup> *<sup>e</sup>* > 0, *N*<sup>∗</sup> *<sup>i</sup>* > 0 and *N*<sup>∗</sup> *<sup>r</sup>* > 0 exist if and only if the parameters of Equations (24)–(27) satisfy

$$R\_2 = \frac{\mu\_2 \theta\_{20} p\_0}{a\_{20} + \gamma\_{20}} < 1,\tag{28}$$

$$R\_3 = \frac{\mu\_3 \theta\_{30} q\_0}{a\_{30} + \gamma\_{30}} < 1,\tag{29}$$

$$R\_4 = \frac{\mu\_3 \theta\_{30} (1 - q\_0) \gamma\_{20} + \gamma\_{20} \gamma\_{30} (\mu\_4 \theta\_{40} + \rho\_0) (a\_{40} + \rho\_0)^{-1}}{(a\_{20} + \gamma\_{20} - \mu\_2 \theta\_{20})(a\_{30} + \gamma\_{30} - \mu\_3 \theta\_{30} q\_0)} < 1,\tag{30}$$

where *R*2, *R*3, and *R*<sup>4</sup> are dimensionless indexes for the existence of the endemic equilibrium point *H*∗. It can be known that *R*<sup>2</sup> < 1 presents that the death rate and conversion rate of exposed population outweigh their reproductive rate. The density of exposes cannot increase indefinitely because of the balance of death rate, conversion rate, and reproductive rate. From the biological point of view, higher death rate and conversion rate are the consequence of infectious disease. *R*<sup>3</sup> < 1 denotes that reproductive rate of infected population is less than their death rate and conversion rate. From the biological point of view, higher death rate and conversion rate in both cases are the fundamental results of infectious diseases. Let *M*<sup>0</sup> = (*α*<sup>20</sup> + *γ*<sup>20</sup> − *μ*2*θ*20)(*α*<sup>30</sup> + *γ*<sup>30</sup> − *μ*3*θ*30*q*0) − *μ*3*θ*30(1 − *q*0)*γ*<sup>20</sup> − *γ*20*γ*30(*μ*4*θ*<sup>40</sup> + *ρ*0)(*α*<sup>40</sup> + *ρ*0)−1. The inequality *R*<sup>4</sup> < 1 is equivalent to *M*<sup>0</sup> > 0. Then, the endemic equilibrium point *H*∗ = (*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) exists and can be explicitly expressed as

$$\begin{split} N\_{\varepsilon}^{\*} &= \frac{(\mu\_{20} + \gamma\_{20} - \mu\_{2}\theta\_{20}p\_{0})(\alpha\_{30} + \gamma\_{30} - \mu\_{3}\theta\_{30}q\_{0})}{\gamma\_{20}\gamma\_{10}\theta\_{0}}, N\_{\varepsilon}^{\*} = \frac{(\mu\_{30} + \gamma\_{30} - \mu\_{3}\theta\_{30}q\_{0})}{\gamma\_{20}}N\_{i}^{\*},\\ N\_{i}^{\*} &= \frac{(\mu\_{1}\theta\_{10} - \mu\_{10})\gamma\_{20}}{M\_{0}}N\_{s}^{\*}, N\_{r}^{\*} = \frac{\gamma\_{30}}{(\mu\_{40} + \rho\_{0})}N\_{i}^{\*}.\end{split}$$

First, the local stability of the disease-free equilibrium point *H*<sup>0</sup> is analyzed.

**Theorem 3.** *For all τ* ≥ 0*, the disease-free equilibrium point H*<sup>0</sup> = (0*,* 0*,* 0*,* 0) *is locally asymptotically stable if R*<sup>1</sup> < 1*, R*<sup>2</sup> < 1 *and R*<sup>3</sup> < 1*.*

**Proof.** The characteristic equation of the system in Equation (20) for the disease-free equilibrium point *H*<sup>0</sup> is

$$(\lambda + \mathfrak{a}\_{10} - \mu\_1 \theta\_{10})(\lambda + \mathfrak{a}\_{20} + \gamma\_{20} - \mu\_2 \theta\_{20} p\_0)(\lambda + \mathfrak{a}\_{30} + \gamma\_{30} - \mu\_3 \theta\_{30} q\_0)(\lambda + \mathfrak{a}\_{40} + \rho\_0) = 0. \tag{31}$$

Four roots *λ*1,2,3,4 of Equation (31) are given:

$$
\lambda\_1 = \mu\_1 \theta\_{10} - \mathfrak{a}\_{10}, \\
\lambda\_2 = -(\mathfrak{a}\_{20} + \gamma\_{20} - \mu\_2 \theta\_{20} p\_0), \\
\lambda\_3 = -(\mathfrak{a}\_{30} + \gamma\_{30} - \mu\_3 \theta\_{30} q\_0), \\
\lambda\_4 = -(\mathfrak{a}\_{40} + \rho\_0).
$$

If *α*<sup>10</sup> > *μ*1*θ*10(*R*<sup>1</sup> < 1), *α*<sup>20</sup> + *γ*<sup>20</sup> > *μ*2*θ*<sup>20</sup> *p*<sup>0</sup> (*R*<sup>2</sup> < 1), and *α*<sup>30</sup> + *γ*<sup>30</sup> > *μ*3*θ*30*q*<sup>0</sup> (*R*<sup>3</sup> < 1), characteristic equation in Equation (31) has four real negative roots. Then, the disease-free equilibrium point *H*<sup>0</sup> = (0, 0, 0, 0) is locally asymptotically stable.

Next, we study the stability of the endemic equilibrium point *H*∗ by linearizing the nonlinear autonomous system and calculating characteristic equations. Accordingly, we can obtain polynomial equations which can analyze the stability of system. After linearization of the system in Equation (20), we get

$$\begin{cases} L\_s'(t) = -(\mathfrak{a}\_{10} - \mu\_1 \theta\_{10}) L\_s(t) - \gamma\_{10} \theta\_0 [L\_s(t) N\_i^\* + N\_s^\* L\_i(t-\tau)] + \mu\_2 \theta\_{20} (1-p\_0) L\_\varepsilon(t) \\ \quad + \mu\_3 \theta\_{30} (1-q\_0) L\_i(t) + (\mu\_4 \theta\_{40} + \rho\_0) L\_\varepsilon(t), \\ L\_\varepsilon'(t) = -(\mathfrak{a}\_{20} + \gamma\_{20} - \mu\_2 \theta\_{20} p\_0) L\_\varepsilon(t) + \gamma\_{10} \theta\_0 [L\_s(t) N\_i^\* + N\_s^\* L\_i(t-\tau)], \\ L\_i'(t) = -(\mathfrak{a}\_{30} + \gamma\_{30} - \mu\_3 \theta\_{30} q\_0) L\_i(t) + \gamma\_{20} L\_\varepsilon(t), \\ L\_\tau'(t) = -(\mathfrak{a}\_{40} + \rho\_0) L\_\tau(t) + \gamma\_{30} L\_i(t). \end{cases} \tag{32}$$

Consider the form of the exponential solution of the linearized system in Equation (32) as

$$L\_s(t) = L\_\varepsilon e^{\lambda \tau},\\ L\_\varepsilon(t) = L\_\varepsilon e^{\lambda \tau},\\ L\_i(t) = L\_i e^{\lambda \tau},\\ L\_r(t) = L\_r e^{\lambda \tau},\tag{33}$$

where *λ* is a parameter. Substituting these solutions into Equation (32), we have

$$\begin{cases} (\lambda + \mathfrak{a}\_{10} - \mu\_1 \theta\_{10} + \gamma\_{10} \theta\_0 N\_i^\*) L\_s + [\gamma\_{10} \theta\_0 N\_s^\* e^{-\lambda \tau} - \mu\_3 \theta\_{30} (1 - q\_0)] L\_i - \mu\_2 \theta\_{20} (1 - p\_0) L\_c \\ - (\mu\_4 \theta\_{40} + \rho\_0) L\_r = 0, \\\ (\lambda + \mathfrak{a}\_{20} + \gamma\_{20} - \mu\_2 \theta\_{20} p\_0) L\_r - \gamma\_{10} \theta\_0 L\_s N\_i^\* - \gamma\_{10} \theta\_0 N\_s^\* L\_i e^{-\lambda \tau} = 0, \\\ (\lambda + \mathfrak{a}\_{30} + \gamma\_{30} - \mu\_3 \theta\_{30} q\_0) L\_i - \gamma\_{20} L\_c = 0, \\\ (\lambda + \mathfrak{a}\_{40} + \rho\_0) \bar{L}\_r - \gamma\_{30} \bar{L}\_i = 0. \end{cases}$$

Denote *M*<sup>1</sup> = *μ*1*θ*<sup>10</sup> − *α*10, *M*<sup>2</sup> = *α*<sup>20</sup> + *γ*<sup>20</sup> − *μ*2*θ*20, *M*<sup>3</sup> = *α*<sup>40</sup> + *ρ*0, *M*<sup>4</sup> = *α*<sup>20</sup> + *γ*<sup>20</sup> − *μ*2*θ*<sup>20</sup> *p*<sup>0</sup> and *M*<sup>5</sup> = *α*<sup>30</sup> + *γ*<sup>30</sup> − *μ*3*θ*30*q*0. Hence *M*<sup>1</sup> > 0, *M*<sup>2</sup> > 0, *M*<sup>3</sup> > 0, *M*<sup>4</sup> > 0 and *M*<sup>5</sup> > 0. The characteristic equation of the system in Equation (32) for the endemic equilibrium point *H*∗=(*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) is given:

$$
\lambda^4 - (\lambda^2 + D\_1 \lambda - D\_2) D\_3 e^{-\lambda \tau} + D\_5 \lambda^3 + D\_6 \lambda^2 + D\_7 \lambda = 0,\tag{34}
$$

where *<sup>D</sup>*<sup>1</sup> <sup>=</sup> *<sup>M</sup>*<sup>3</sup> <sup>−</sup> *<sup>M</sup>*1, *<sup>D</sup>*<sup>2</sup> <sup>=</sup> *<sup>M</sup>*1*M*3, *<sup>D</sup>*<sup>3</sup> <sup>=</sup> *<sup>M</sup>*4*M*5, *<sup>D</sup>*<sup>4</sup> <sup>=</sup> *<sup>M</sup>*<sup>4</sup> <sup>+</sup> *<sup>M</sup>*5, *<sup>D</sup>*<sup>5</sup> <sup>=</sup> *<sup>D</sup>*<sup>1</sup> <sup>+</sup> *<sup>D</sup>*<sup>4</sup> <sup>+</sup> *<sup>M</sup>*1*D*3*M*−<sup>1</sup> <sup>0</sup> , *D*<sup>6</sup> = *<sup>D</sup>*1*D*<sup>4</sup> <sup>−</sup> *<sup>D</sup>*<sup>2</sup> <sup>+</sup> *<sup>D</sup>*<sup>3</sup> + [*D*<sup>4</sup> <sup>−</sup> *<sup>μ</sup>*2*θ*20(<sup>1</sup> <sup>−</sup> *<sup>p</sup>*0) + *<sup>M</sup>*3]*M*1*D*3*M*−<sup>1</sup> <sup>0</sup> , *D*<sup>7</sup> = *D*1*D*<sup>3</sup> − *D*2*D*<sup>4</sup> + [(*M*<sup>2</sup> + *M*5)*M*<sup>3</sup> + *<sup>M</sup>*2*M*<sup>5</sup> <sup>−</sup> *<sup>μ</sup>*3*θ*30(<sup>1</sup> <sup>−</sup> *<sup>q</sup>*0)*γ*20]*M*1*D*3*M*−<sup>1</sup> <sup>0</sup> .

**Theorem 4.** *For the system in Equation (20), if <sup>D</sup>*5(*D*<sup>6</sup> <sup>−</sup> *<sup>D</sup>*3)(*D*<sup>7</sup> <sup>−</sup> *<sup>D</sup>*1*D*3) <sup>−</sup> *<sup>D</sup>*<sup>2</sup> <sup>5</sup>*D*2*D*<sup>3</sup> <sup>−</sup> (*D*<sup>7</sup> <sup>−</sup> *<sup>D</sup>*1*D*3)<sup>2</sup> <sup>&</sup>gt; 0, *the endemic equilibrium point H*∗*=*(*N*∗ *<sup>s</sup> ,N*<sup>∗</sup> *<sup>e</sup> ,N*<sup>∗</sup> *<sup>i</sup> ,N*<sup>∗</sup> *<sup>r</sup>* ) *is locally asymptotically stable when τ* = 0*.*

**Proof.** When *τ* = 0, the characteristic equation (Equation (34)) becomes

$$
\lambda^4 + D\mathfrak{s}\lambda^3 + (D\mathfrak{s} - D\mathfrak{s})\lambda^2 + (D\mathfrak{\gamma} - D\_1D\mathfrak{s})\lambda + D\_2D\mathfrak{s} = 0. \tag{35}
$$

It is evident that *D*<sup>2</sup> > 0, *D*<sup>3</sup> > 0 and *D*<sup>4</sup> > 0 when *N*<sup>∗</sup> *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* satisfy the nonnegativity conditions in Equation (23)–(30). The following equations are also satisfied:

$$\begin{aligned} D\_3 M\_0^{-1} - 1 &= [\mu\_3 \theta\_{30} (1 - q\_0) \gamma\_{20} - \gamma\_{20} \gamma\_{30} (\mu\_4 \theta\_{40} + \rho\_0) M\_3^{-1}] M\_0^{-1} > 0, \\ D\_5 &= D\_4 + M\_1 (D\_3 M\_0^{-1} - 1) + M\_3 > 0, \\ D\_6 - D\_3 &= D\_2 (D\_3 M\_0^{-1} - 1) + M\_1 M\_5 D\_3 M\_0^{-1} + M\_4 [\mu\_3 \theta\_{30} (1 - q\_0) \gamma\_{20} + \gamma\_{20} \gamma\_{30} (\mu\_4 \theta\_{40} + \rho\_0) M\_3^{-1}] \\ &> 0, \\ D\_7 - D\_1 D\_3 &= M\_4 [\mu\_3 \theta\_{30} (1 - q\_0) \gamma\_{20} + \gamma\_{20} \gamma\_{30} (\mu\_4 \theta\_{40} + \rho\_0) M\_3^{-1}] D\_2 M\_0^{-1} + M\_5 D\_2 (D\_3 M\_0^{-1} - 1) \\ &+ M\_1 D\_3 + \gamma\_{20} \gamma\_{30} (\mu\_4 \theta\_{40} + \rho\_0) M\_3^{-1} M\_1 D\_3 M\_0^{-1} > 0, \\ D\_2 D\_3 &= M\_1 M\_3 M\_4 M\_5 > 0. \end{aligned}$$

Suppose that

$$D\_5(D\_6 - D\_3)(D\_7 - D\_1 D\_3) - D\_5^2 D\_2 D\_3 - (D\_7 - D\_1 D\_3)^2 > 0,\tag{36}$$

Then, according to the Routh–Hurwitz criterion, the roots of the characteristic equations of system are all negative real. Thus, the endemic equilibrium point of the autonomous system in Equation (20) is asymptotically stable with *τ* = 0.

Next, the characteristic equation (Equation (34)) with *τ* > 0 is presented for discussion.

**Theorem 5.** *For the system in Equation (20), if D*<sup>2</sup> <sup>6</sup> <sup>−</sup> <sup>2</sup>*D*5*D*<sup>7</sup> <sup>−</sup> *<sup>D</sup>*<sup>2</sup> <sup>3</sup> > <sup>0</sup> *and <sup>D</sup>*<sup>2</sup> <sup>7</sup> <sup>−</sup> <sup>2</sup>*D*2*D*<sup>2</sup> <sup>3</sup> <sup>−</sup> *<sup>D</sup>*<sup>2</sup> 1*D*<sup>2</sup> <sup>3</sup> > 0*, the results can be obtained:*

(i) When 0 < *τ* < *τ*0, the endemic equilibrium point *H*∗=(*N*<sup>∗</sup> *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) is locally asymptotically stable.

(ii) When *τ* > *τ*0, the endemic equilibrium point *H*∗=(*N*<sup>∗</sup> *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) is unstable.

(iii) A Hopf bifurcation occurs at *H*∗ = (*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) when *τ* = *τ<sup>k</sup>* (*k* = 0, 1, 2, ···).

**Proof.** When *τ* > 0, suppose that Equation (34) has a purely imaginary root *λ* = *iω*(*ω* > 0). The characteristic equation (Equation (34)) converses to the following form:

$$-\omega^4 - (-\omega^2 + i\mathcal{D}\_1\omega - \mathcal{D}\_2)\mathcal{D}\_3(\cos(\omega\tau) - i\sin(\omega\tau)) - i\mathcal{D}\_5\omega^3 - \mathcal{D}\_6\omega^2 + i\mathcal{D}\_7\omega = 0.$$

Separating the real and imaginary parts yields two corresponding equations:

$$
\omega D\_1 D\_3 \sin(\omega \tau) - (\omega^2 D\_3 + D\_2 D\_3) \cos(\omega \tau) = \omega^4 - D\_6 \omega^2,\tag{37}
$$

$$(\omega^2 D \mathfrak{z} + D \mathfrak{z} D \mathfrak{z}) \sin(\omega \tau \mathfrak{x}) + \omega \mathfrak{z} D\_1 D \mathfrak{z} \cos(\omega \tau \mathfrak{x}) = D \mathfrak{z} \omega - D \mathfrak{z} \omega^3. \tag{38}$$

Thus,

$$\sin\omega\tau = \frac{(\omega^2 D\_3 + D\_2 D\_3)(D\_7 \omega - D\_5 \omega^3) + \omega D\_1 D\_3 (\omega^4 - D\_6 \omega^2)}{(\omega D\_1 D\_3)^2 + (\omega^2 D\_3 + D\_2 D\_3)^2},\tag{39}$$

$$\cos\omega\tau = \frac{\omega D\_1 D\_3 (D\_7 \omega - D\_5 \omega^3) - (\omega^4 - D\_6 \omega^2)(\omega^2 D\_3 + D\_2 D\_3)}{(\omega D\_1 D\_3)^2 + (\omega^2 D\_3 + D\_2 D\_3)^2}. \tag{40}$$

Let *x* = *ω*2. Squaring Equations (37) and (38) and adding them together, we can get

$$\mathbf{x}^4 + (D\_5^2 - 2D\_6)\mathbf{x}^3 + (D\_6^2 - 2D\_5D\_7 - D\_3^2)\mathbf{x}^2 + (D\_7^2 - 2D\_2D\_3^2 - D\_1^2D\_3^2)\mathbf{x} - D\_2^2D\_3^2 = 0. \tag{41}$$

Let *l*(*x*) = *x*<sup>4</sup> + (*D*<sup>2</sup> <sup>5</sup> <sup>−</sup> <sup>2</sup>*D*6)*x*<sup>3</sup> + (*D*<sup>2</sup> <sup>6</sup> <sup>−</sup> <sup>2</sup>*D*5*D*<sup>7</sup> <sup>−</sup> *<sup>D</sup>*<sup>2</sup> 3)*x*<sup>2</sup> + (*D*<sup>2</sup> <sup>7</sup> <sup>−</sup> <sup>2</sup>*D*2*D*<sup>2</sup> <sup>3</sup> <sup>−</sup> *<sup>D</sup>*<sup>2</sup> 1*D*<sup>2</sup> <sup>3</sup>)*<sup>x</sup>* <sup>−</sup> *<sup>D</sup>*<sup>2</sup> 2*D*<sup>2</sup> 3. According to the previous restrictions in Equation (23)–(30), we can also get

$$D\_5^2 - 2D\_6 = M\_3^2 + M\_4^2 + M\_5^2 + M\_1^2(D\_3M\_0^{-1} - 1)^2 + 2\mu\_2\theta\_{20}(1 - p\_0)M\_1D\_3M\_0^{-1} > 0. \tag{42}$$

Suppose that

$$D\_6^2 - 2D\_5D\_7 - D\_3^2 > 0,\tag{43}$$

$$D\_7^2 - 2D\_2D\_3^2 - D\_1^2D\_3^2 > 0.\tag{44}$$

In accordance with Descartes' rule of signs, Equation (41) has only one changed sign. Thus, it only has one positive root. Let *x*∗ be the small unique positive root and it always exists. The unknown parameter *ω* of Equations (37) and (38) is defined as ±*iω*<sup>0</sup> = ±*i* <sup>√</sup>*x*0. Combining with Equations (37) and (38), the form of time delay *τ<sup>k</sup>* is gained:

$$\pi\_k = \frac{1}{\omega\_0} \arcsin(\frac{(\omega\_0^2 D\_3 + D\_2 D\_3)(D\_7 \omega\_0 - D\_5 \omega\_0^3) + \omega\_0 D\_1 D\_3 (\omega\_0^4 - D\_6 \omega\_0^2)}{(\omega\_0 D\_1 D\_3)^2 + (\omega\_0^2 D\_3 + D\_2 D\_3)^2} + 2k\pi), k = 0, 1, 2, \cdots \tag{45}$$

Therefore, when *τ* ∈(0, *τ*0), all roots of Equation (34) have strictly negative real parts. The endemic equilibrium point *H*∗=(*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) of the system in Equation (20) is locally asymptotically stable. When *τ* = *τ*0, the roots of Equation (34) have strictly negative real parts except for ±*iω*0. When *τ* > *τ*0, the endemic equilibrium point *H*∗ = (*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) of the system in Equation (20) is unstable. Then, differentiating both sides of Equation (34) with respect to *τ*, we obtain

$$(\frac{d\lambda}{d\tau})^{-1} = \frac{(2\lambda + D\_1)D\_3 - (4\lambda^3 + 3D\_5\lambda^2 + 2D\_6\lambda + D\_7)e^{\lambda\tau}}{(\lambda^3 + D\_1\lambda^2 - D\_2\lambda)D\_3} - \frac{\tau}{\lambda}. \tag{46}$$

Further, one leads to

$$\text{Re}\left\{ (\frac{d\lambda}{d\tau})^{-1} \right\}\_{\lambda = i\omega\_0} = \frac{Q\_1 + Q\_2 - Q\_3}{\omega\_0^2 [(\omega\_0 D\_1 D\_3)^2 + (\omega\_0^2 D\_3 + D\_2 D\_3)^2]} \tag{47}$$

where

$$\begin{aligned} Q\_1 &= \omega\_0^2 D\_1 D\_3 (4\omega\_0^3 - 2D\_6 \omega\_0) \sin \omega\_0 \tau\_0 + \omega\_0^2 D\_1 D\_3 (D\_7 - 3D\_5 \omega\_0^2) \cos \omega\_0 \tau\_0, \\ Q\_2 &= \omega\_0 (\omega\_0^2 D\_3 + D\_2 D\_3) (D\_7 - 3D\_5 \omega\_0^2) \sin \omega\_0 \tau\_0 - \omega\_0 (\omega\_0^2 D\_3 + D\_2 D\_3) (4\omega\_0^3 - 2D\_6 \omega\_0) \cos \omega\_0 \tau\_0, \\ Q\_3 &= 2\omega\_0^2 D\_3 (\omega\_0^2 D\_3 + D\_2 D\_3) + \left(\omega \eta D\_1 D\_3\right)^2. \end{aligned}$$

Since conditions in Equation (42)–(44) hold, then

$$\operatorname{sign}\left\{\operatorname{Re}(\frac{d\lambda}{d\tau})\right\}\_{\tau=\tau\_0} = \operatorname{sign}\left\{\operatorname{Re}(\frac{d\lambda}{d\tau})^{-1}\right\}\_{\lambda=i\omega\_0} = \operatorname{sign}\left\{\frac{l'(\mathbf{x}\_0)}{(\omega\_0 D\_1 D\_3)^2 + (\omega\_0^2 D\_3 + D\_2 D\_3)^2}\right\}\_{\mathbf{x}\_0 = \omega\_0^2} > 0,\tag{48}$$

where *l* (*x*0) = 4*x*<sup>3</sup> <sup>0</sup> + <sup>3</sup>(*D*<sup>2</sup> <sup>5</sup> <sup>−</sup> <sup>2</sup>*D*6)*x*<sup>2</sup> <sup>0</sup> + 2(*D*<sup>2</sup> <sup>6</sup> <sup>−</sup> <sup>2</sup>*D*5*D*<sup>7</sup> <sup>−</sup> *<sup>D</sup>*<sup>2</sup> <sup>3</sup>)*x*<sup>0</sup> + (*D*<sup>2</sup> <sup>7</sup> <sup>−</sup> <sup>2</sup>*D*2*D*<sup>2</sup> <sup>3</sup> <sup>−</sup> *<sup>D</sup>*<sup>2</sup> 1*D*<sup>2</sup> <sup>3</sup>) > 0. Hence, based on the properties of the Hopf bifurcation discussed in [32], the transversal condition holds and a Hopf bifurcation occurs at *τ* = *τ*0. The proof is complete.

#### **5. Simulation**

To affirm the stability analysis above, we numerically simulated the disease-free equilibrium point *H*<sup>0</sup> and the endemic equilibrium point *H*∗ of the system in Equation (20). We are more concerned with the indicators and conditions for outbreaks of "local" population or local asymptotic stability of equilibrium points. All programs were developed using the software Matlab R2016a (see Program S1).

First, we considered the stability of the disease-free equilibrium point *H*<sup>0</sup> = (0, 0, 0, 0). Let *αi*<sup>0</sup> = 0.13, *i* = 1, 2, 3, 4; *μ*<sup>1</sup> = 0.3; *μ*<sup>2</sup> = 0.5; *μ*<sup>3</sup> = 0.5; *μ*<sup>4</sup> = 0.3; *θ*<sup>10</sup> = 0.4; *θ*<sup>20</sup> = 0.1; *θ*<sup>30</sup> = 0.1; *θ*<sup>40</sup> = 0.2; *γ*<sup>10</sup> = 0.1; *γ*<sup>20</sup> = 0.15; *γ*<sup>30</sup> = 0.1; *β*<sup>0</sup> = 0.0023; *p*<sup>0</sup> = 0.1; *q*<sup>0</sup> = 0.15; *ρ*<sup>0</sup> = 0.02; and the initial value was (*Ns*(0), *Ne*(0), *Ni*(0), *Nr*(0)) = (1500, 1000, 500, 0). The parameters of the system in Equation (20) satisfy conditions *R*<sup>1</sup> < 1, *R*<sup>2</sup> < 1, and *R*<sup>3</sup> < 1 and the conditions in Theorem 3 are satisfied in Figure 1. The disease-free equilibrium point *H*<sup>0</sup> of the system in Equation (20) is locally asymptotically stable.

**Figure 1.** Local asymptotic stability of the disease-free equilibrium point *H*<sup>0</sup> = (0, 0, 0, 0).

Next, we considered the endemic equilibrium point *H*∗ = (*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) with *τ* = 0. Let *αi*<sup>0</sup> = 0.13, *i* = 1, 2, 3, 4; *μ*<sup>1</sup> = 0.6; *μ*<sup>2</sup> = 0.5; *μ*<sup>3</sup> = 0.5; *μ*<sup>4</sup> = 0.3; *θ*<sup>10</sup> = 0.4; *θ*<sup>20</sup> = 0.1; *θ*<sup>30</sup> = 0.1; *θ*<sup>40</sup> = 0.2; *γ*<sup>10</sup> = 0.1; *γ*<sup>20</sup> = 0.15; *γ*<sup>30</sup> = 0.1; *β*<sup>0</sup> = 0.0023; *p*<sup>0</sup> = 0.1; *q*<sup>0</sup> = 0.15; *ρ*<sup>0</sup> = 0.02; and the initial value was (*Ns*(0), *Ne*(0), *Ni*(0), *Nr*(0)) = (1500, 1000, 500, 0). The parameters of the system in Equation (20) satisfy the condition in Equation (36). The result is presented in Figure 2, which comes from Theorem 4; the system in Equation (20) has a unique endemic equilibrium point *H*∗ = (*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* )≈ (1723, 1108, 797, 531) and the system in Equation (20) is locally asymptotically stable.

**Figure 2.** Local asymptotic stability of the endemic equilibrium point *H*∗ = (*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) when *τ* = 0.

Similarly, let *αi*<sup>0</sup> = 0.13, *i* = 1, 2, 3, 4; *μ*<sup>1</sup> = 0.6; *μ*<sup>2</sup> = 0.5; *μ*<sup>3</sup> = 0.5; *μ*<sup>4</sup> = 0.3; *θ*<sup>10</sup> = 0.4; *θ*<sup>20</sup> = 0.1; *θ*<sup>30</sup> = 0.1; *θ*<sup>40</sup> = 0.2; *γ*<sup>10</sup> = 0.1; *γ*<sup>20</sup> = 0.16; *γ*<sup>30</sup> = 0.1; *β*<sup>0</sup> = 0.0023; *p*<sup>0</sup> = 0.1; *q*<sup>0</sup> = 0.15; *ρ*<sup>0</sup> = 0.02; and the initial value be (*Ns*(0), *Ne*(0), *Ni*(0), *Nr*(0))=(1500, 1000, 500, 0). In other words, *Ns*(0) < *N*<sup>∗</sup> *s* , *Ne*(0) < *N*<sup>∗</sup> *<sup>e</sup>* , *Ni*(0) < *N*<sup>∗</sup> *<sup>i</sup>* , and *Nr*(0) < *N*<sup>∗</sup> *<sup>r</sup>* . The parameters of the system in Equation (20) satisfy the conditions in Equations (43) and (44). In this case, we have the roots of Equation (41) are *x*<sup>1</sup> = −0.1398, *x*2=0.0073, and *x*3,4 = −0.213 ± 0.0248*i*. Thus, we get *x*<sup>0</sup> = 0.0073, *ω*<sup>0</sup> = 0.0855, *τ*<sup>0</sup> = 5.3526, and *τ*<sup>1</sup> = 73.4875. The population dynamic behaviors of the endemic equilibrium point *H*<sup>∗</sup> are shown in Figure 3 (*τ* = 4.3526) and Figure 4 (*τ* = 6.3526), respectively. It can been shown that the endemic equilibrium point *H*<sup>∗</sup> spends longer time to become locally asymptotically stable only when 0 < *τ* < *τ*<sup>0</sup> (see Figure 3). When *τ* > *τ*0, it is quite clear that larger values of the time delay causes periodic oscillations. The equilibrium point loses its stability. Thus, the solution of the system in Equation (20) is unstable with the increase of time t (see Figure 4).

**Figure 3.** Local asymptotic stability of the endemic equilibrium point *H*∗ = (*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) when *τ* = 4.3526 < *τ*0. (**a**) Time response diagram of *Ns*, *Ne*, *Ni* and *Nr*; (**b**) Phase diagram of *Ns*, *Ne* and *Ni*.

**Figure 4.** Instability of the endemic equilibrium point *H*∗ = (*N*∗ *<sup>s</sup>* , *N*<sup>∗</sup> *<sup>e</sup>* , *N*<sup>∗</sup> *<sup>i</sup>* , *N*<sup>∗</sup> *<sup>r</sup>* ) when *τ* = 6.3526 > *τ*0. (**a**) Time response diagram of *Ns*; (**b**) Time response diagram of *Ne*; (**c**) Time response diagram of *Ni*; (**d**) Time response diagram of *Nr*; (**e**) Phase diagram of *Ns*, *Ne* and *Ni*.
