Recurrent epidemic waves in a delayed epidemic model with quarantine

ABSTRACT In this paper, we are concerned with an epidemic model with quarantine and distributed time delay. We define the basic reproduction number and show that if , then the disease-free equilibrium is globally asymptotically stable, whereas if , then it is unstable and there exists a unique endemic equilibrium. We obtain sufficient conditions for a Hopf bifurcation that induces a nontrivial periodic solution which represents recurrent epidemic waves. By numerical simulations, we illustrate stability and instability parameter regions. Our results suggest that the quarantine and time delay play important roles in the occurrence of recurrent epidemic waves.


Introduction
In the Coronavirus Disease Outbreak 2019 (COVID- 19), mathematical models have been widely applied for epidemiological considerations (see, e.g. [7,8]). In COVID-19, recurrent epidemic waves have been observed in many countries [20]. One of the purposes of this paper is to discuss the mechanism of such recurrent epidemic waves from the viewpoint of mathematical modelling. In COVID-19, strict non-pharmaceutical intervention policies such as lockdown were taken to suppress the epidemic in many countries. Effects of such policies have been estimated in many studies, and they are often concluded to be positive in terms of disease suppression (see, e.g. [9,16]). However, the rebound of epidemic after the lifting of strict interventions was observed in many countries [20]. We can guess that people's behaviour change induced by such strict interventions is one of the factors of recurrent epidemic waves. For this perspective, in this paper, we will construct an epidemic model with quarantine and study the existence of periodic solutions which represent the recurrent epidemic waves.
Periodicity in infectious disease outbreaks has been studied for decades from the viewpoint of mathematical modelling [13]. Periodic coefficients, time delay, nonlinear incidence, variable population size and age structure [6] are known to be the factors of periodic solutions in epidemic models [13]. In this paper, we are concerned with periodic solutions caused by time delay in an epidemic model with a quarantined class.
Periodic solutions of delayed epidemic models have been studied in many papers (see, e.g. [1][2][3]15,17,19,[21][22][23]). In many papers, models with constant time delay were considered, in which time delay itself plays the role of a bifurcation parameter for a Hopf bifurcation. However, distributed time delay (see, e.g. [4]) might be more realistic because it enables us to consider the effect of the history of solutions in a fixed time interval. In this paper, we construct an epidemic model with distributed time delay, and obtain sufficient conditions for a Hopf bifurcation which induces nontrivial periodic solutions.
One of the most important concepts in mathematical epidemiology is the basic reproduction number R 0 [14,Chapter 9], which is the expected number of secondary cases produced by a typical infectious individual in a full susceptible population [5]. Intuitively speaking, if R 0 < 1, then there is no outbreak, whereas if R 0 > 1, then an outbreak occurs. In epidemic models, R 0 < 1 (resp., R 0 > 1) often implies the global stability of the disease-free (resp., endemic) equilibrium [14,Section 5.5.2]. However, in some epidemic models, global stability of the disease-free (resp., endemic) equilibrium for R 0 < 1 (resp., R 0 > 1) is not true by virtue of characteristic phenomena such as a backward bifurcation at R 0 = 1 [11] and a Hopf bifurcation [10]. Therefore, investigation of the relation between R 0 and the stability of each equilibrium in epidemic models is important nontrivial work. Such investigation would help us to understand the role of R 0 in controlling the epidemic dynamics. In this paper, we investigate the relation between R 0 and the stability of each equilibrium of our model.
The organization of this paper is as follows: in Section 2, we construct our main model with quarantine and distributed time delay. We show the well-posedness of the problem and define the basic reproduction number R 0 . In Section 3, we show that if R 0 ≤ 1, then the disease-free equilibrium of our model is globally asymptotically stable, whereas if R 0 > 1, then it is unstable. Moreover, we obtain sufficient conditions for the local asymptotic stability of the endemic equilibrium when R 0 > 1. These conditions correspond to the situation where the quarantine is less restrictive. In Section 4, we study the characteristic equation of the endemic equilibrium for R 0 > 1 and obtain sufficient conditions for a Hopf bifurcation. In Section 5, we give numerical examples that exhibit periodic solutions induced by a Hopf bifurcation. Our results suggest that the quarantine and time delay play important roles in the occurrence of recurrent epidemic waves. Finally, Section 6 is devoted to a discussion.

Model and preliminaries
Let S(t), Q(t), I(t) and R(t) be the susceptible, quarantined, infectious and removed populations at time t, respectively. Let b be the birth rate, μ be the mortality rate, β be the infection rate and γ be the removal rate. Let g be the quarantine rate at which susceptible individuals go into the quarantined class, and h be the lifting rate at which quarantined individuals go back to the susceptible class. We assume that g and h depend on the history of the infectious population, that is, g = g(I) and h = h(I), which are given as follows: where q and α represent the sensitivity of the quarantine and lifting rates to the history of the infectious population, respectively. δ is the lifting rate when and f (τ ) represents the intensity of how the infectious population at τ time ago affects the quarantine and lifting rates at the current time. Here, we can regard U(t) := ∞ 0 f (τ )I(t − τ ) dτ as the amount of valid information on the history of the infectious population at time t. We assume that the quarantine rate g and the lifting rate h are proportional and inversely proportional to U, respectively. The main model in this paper is given by the following system of delay differential equations (see also Figure 1): We make the following assumptions on each parameter: (A1) b, μ, β, γ , q, δ and α are strictly positive. (A2) 0 < σ < 1.
As equations of S, Q and I in system (1) are independent of R, it is sufficient to consider the following reduced system: Let E be the set of all bounded and continuous functions from (−∞, 0] to R 3 , equipped with the norm where · ∞ denotes the maximum norm in R 3 . That is, r ∞ = max(|r 1 |, |r 2 |, |r 3 |) for all r = (r 1 , r 2 , r 3 ) T ∈ R 3 , where T denotes the transpose operation. We define the positive cone of E by E + := E ∩ C((−∞, 0], R 3 + ), and the solution space for system (2) by where N = b/μ. For any t ∈ R, we write u(t) := (S(t), Q(t), I(t)) T . Moreover, for any t ≥ 0 and ρ ∈ (−∞, 0], we write u t (ρ) := u(t + ρ). The system (2) can then be rewritten as Let φ := (φ 1 , φ 2 , φ 3 ) T be the initial condition for system (2), that is, for all ρ ≤ 0. We now prove the following proposition.

Proposition 2.1:
is positively invariant for system (2). That is, if φ ∈ , then u t ∈ for all t > 0.
Proof: We first show that F := (F 1 , F 2 , F 3 ) T is Lipschitz continuous in . Indeed, for any ϕ, η ∈ , we have In a similar manner, we can show that there exist positive constants K 2 and K 3 , which are independent of ϕ and η, such that 3. This implies that F is Lipschitz continuous in . By the Lipschitz continuity of F, the existence and uniqueness of a local solution u for φ = (φ 1 , φ 2 , φ 3 ) T ∈ follows from [12, Theorem 2.3 in Chapter 2]. The positivity of I is obvious as Suppose that there exists a t 1 > 0 such that S(t) ≥ 0 and Q(t) ≥ 0 for all t ∈ (0, t 1 ) and S(t 1 ) = 0. From the first equation of (2), we have S (t 1 ) = b + h(I)(t 1 )Q(t 1 ) ≥ b > 0, which implies the positivity of S. In a similar manner, we can show the positivity of Q.
Hence, we can define the basic reproduction number R 0 [5] for system (2) by implies the exponential decay (resp., growth) of the infectious population near the disease-free equilibrium E 0 . More precisely, we can show that E 0 is globally asymptotically stable if R 0 ≤ 1, whereas E 0 is unstable if R 0 > 1 (see Propositions 3.1 and 3.2 in the next section). This implies that the backward bifurcation at R 0 = 1 never occurs in our model and reducing R 0 less than 1 is a clear goal to control the epidemic. The next proposition states that R 0 also determines the existence of a unique endemic equilibrium E * in .

Proposition 2.2:
If R 0 > 1, then system (2) has a unique endemic equilibrium E * in . If R 0 ≤ 1, then system (2) does not have any endemic equilibrium E * in . Proof: One can see from (2) and (A3) that E * : (S, Q, I) = (S * , Q * , I * ) is a solution to the following system of algebraic equations: (3) By the third equation of (3), we have Moreover, by adding the three equations of (3), we have Substituting (4) into (5) gives On the other hand, substituting (4) into the second equation of (3) gives and hence, Substituting (7) into (6), we have Rearranging (8), we arrive at the cubic equation Suppose that R 0 > 1. As c 3 > 0 and G(0) = c 0 < 0, G(x) = 0 has at least one positive root I * > 0. One can easily see from (8) that I * < b/(μ + γ ) < N. Moreover, Q * > 0 obviously follows from (7). One can then see from (3) that Consequently, E * exists in . The uniqueness follows because the right-hand side of (8) is monotone decreasing with respect to I * .

Stability of equilibria
Let X = S −S, W = Q −Q and Y = I −Ī be the perturbations from the disease-free or endemic equilibrium (S,Q,Ī) T ∈ . We then have Hence, the system (2) is linearized around (S,Q,Ī) T as follows: The following proposition states that the local stability of the disease-free equilibrium E 0 is determined by the basic reproduction number R 0 . Proposition 3.1: If R 0 < 1, then the disease-free equilibrium E 0 : (S, Q, I) = (N, 0, 0) of system (2) is locally asymptotically stable. If R 0 > 1, then E 0 is unstable.
More strongly, we can show the global asymptotic stability of the disease-free equilibrium E 0 for R 0 ≤ 1 by using a construction method of Lyapunov functions in [18]. (2) is globally asymptotically stable in .

Proposition 3.2:
Proof: Inspired by Lyapunov functions in [18], we define The derivative of V along the solution trajectory of (2) is then calculated as We then see that , we see that the disease-free equilibrium E 0 is globally asymptotically stable in . This completes the proof.
Proposition 3.2 implies that to make R 0 less than 1 is a clear goal to curb the epidemic in our model. In what follows, we focus on the case where R 0 > 1. As shown in Proposition 2.2, there exists a unique endemic equilibrium E * of system (2). By (4) and (9), the characteristic equation for E * is given by We next consider the following two cases: (i) q is small; (ii) δ is large. These cases represent the situation where the quarantine is less restrictive. The next proposition states that E * is asymptotically stable in these cases.

Proposition 3.3:
Suppose that R 0 > 1. The endemic equilibrium E * of system (2) is locally asymptotically stable if either (i) q is sufficiently small, or (ii) δ is sufficiently large.
Proof: (i) By the continuity, it suffices to consider the case where q → +0. Since 0 < S * , I * < N, qS * , qI * → +0 as q → +0. Moreover, by (4), (7) and (8), one can see that Consequently, It is clear that the cubic equation in the right-hand side only has roots with negative real parts. This implies that E * is asymptotically stable if q is sufficiently small.
(ii) By the continuity, it suffices to consider the case where δ → +∞. Similarly, as in (i), by (4), (7) and (8), one can see that Moreover, by (7), one can see that Multiplying 1/δ by P(λ) gives As δ → +∞, this converges to It is obvious that the quadratic equation in the right-hand side only has roots with negative real parts. This implies that E * is asymptotically stable if δ is sufficiently large. The proof is complete.
Proposition 3.3 suggests that if the quarantine is not so restrictive, then no recurrent epidemic wave can occur because the endemic equilibrium is asymptotically stable.
The characteristic equation P(λ) = 0 for E * can then be rearranged as follows: where a 2 = [q + (1 + σ )β]I * + 2μ +δ, We now assume that f is given by the following truncated uniform distribution: where T ≥ 0 denotes the time delay (e.g. incubation period) until infectious individuals are reported and L > 0 denotes the length of the period in which the information of the past reported individuals affect the current people's behaviour. The characteristic Equation (11) can then be rewritten as Note that λ = 0 because P(0) = a 0 + b 0 > 0. The next proposition states that the endemic equilibrium E * is asymptotically stable if T and L are sufficiently small, i.e. there is less time delay effect.

Proposition 4.1:
Suppose that R 0 > 1 and f is given by (12). If both T and L are sufficiently small, then the endemic equilibrium E * of system (2) is locally asymptotically stable.
Proof: Letting T, L → +0 in (13), we have Note that all coefficients in the left-hand side are strictly positive. Hence, by the Routh-Hurwitz criterion, it suffices to show that a 2 (a 1 + b 1 ) − (a 0 + b 0 ) > 0. In fact, it is easy to see that We are now in a position to show that a 2 a 1 > a 0 . By (4) and (7), we have (15) Using (15), we have By (14) and (16), we obtain a 2 (a 1 + b1) − (a 0 + b 0 ) > 0. This completes the proof.
Proposition 4.1 suggests that there can be no recurrent epidemic wave if there is less time delay in the information on the infectious population. To study the Hopf bifurcation of E * , we seek the existence of a conjugate pair of purely imaginary roots λ = ±iω (ω > 0) of the characteristic Equation (13). Rearranging (13) and taking absolute values, we have Substituting λ = ±iω and taking squares of both sides, we have Thus, we obtain Letting x = ω 2 , this equation can be rewritten as One can then see that if h(x) = 0 has a positive root x * > 0, then λ = ±i √ x * is a conjugate pair of purely imaginary roots of the characteristic Equation (13). Indeed, the following proposition states that the direction of the crossing of the imaginary axis is determined by the sign of h (x * ).

Proposition 4.2:
Suppose that R 0 > 1 and f is given by (12). If there exists an x * > 0 such that h(x * ) = 0, then the characteristic Equation (13) has a conjugate pair λ = ±i √ x * of purely imaginary roots. Moreover, if h (x * ) > 0 (resp., h (x * ) < 0), then the pair λ = ±i √ x * crosses the imaginary axis from left to right (resp. right to left) as T increases.
Proof: The first part follows from the above argument. To show the last part, we rewrite the characteristic Equation (13) as Let λ = λ(T) and differentiating (18) with respect to T gives Rearranging gives As the sign of Re λ is equal to that of Re(1/λ ), it suffices to check the sign of Re(1/λ ) at λ = iω to judge the direction of the crossing of the imaginary axis. That is, if Re (1/λ ) > 0 (resp., < 0) at λ = iω, then the conjugate pair λ = ±iω crosses the imaginary axis from left to right (resp., right to left). At λ = iω, we have and thus, by (17), .
As x * = ω 2 , we see that the sign of Re (1/λ ) is equal to that of h (x * ). This completes the proof.
Under the assumptions in Proposition 4.2, the Hopf bifurcation of the endemic equilibrium E * can occur regarding T as a bifurcation parameter. Let T * be a critical value such that λ(T * ) = iω = i √ x * satisfies the characteristic Equation (13). That is, By Euler's formula, we have e −iωT * = cos ωT * − i sin ωT * , and thus, More precisely, we have , .
Hence, we can compute T * as In the next section, using these results, we will provide numerical examples in which the Hopf bifurcation arises, that is, the recurrent epidemic waves occur.

Numerical experiments
We fix the following parameters for system (2) with distribution (12): Although these parameters were not estimated from any real data, the epidemiological justification of the choice is as follows: the unit time is one year, and • The total population N = b/μ is normalized to 1.
• The basic reproduction number R 0 is 2.5.
• The length of the validity period of information is L = 1/4 years = 3 months.
In Figure 3, T = 3.8 ∈ T s . In this case, the infectious population I converges to the steady state I * ≈ 3.121 × 10 −4 as time evolves, which implies that E * is asymptotically stable.
In Figure 4, T = 4.2 ∈ T u . In this case, the infectious population I converges to a nontrivial periodic solution as time evolves, which implies that E * is unstable.
We see from Figure 5 that the bifurcation at T = T + n , n = 0, 1, 2 is supercritical, whereas the bifurcation at T = T − n , n = 0, 1, 2 is subcritical. One might think that T + 0 ≈ 3.9919 years is unrealistic for the time delay of information that causes the recurrent epidemic waves. In fact, we can obtain more plausible values of the parameter interval (T + 0 , T − 0 ) by changing the sensitivity q and α (see Table 1). We see from Table 1 that the interval (T + 0 , T − 0 ) is enlarged as the sensitivity q and α increase. For example, if q = α = 10,000, then T + 0 ≈ 0.0048 years ≈ 1.752 days.   For fixed α = 1, the stability and instability regions depicted by curves T + n and T − n , n = 0, 1, 2 with varied q are shown in Figure 6(a).
One can see from this figure that the endemic equilibrium E * is asymptotically stable when q is small, and it is destabilized as q becomes large. This result is consistent with  Figure 6. Stability and instability regions depicted by curves T + n and T − n , n = 0, 1, 2 for parameters (19) and α = 1. (a) 0 ≤ q ≤ 500 (δ = 12); (b) 0 ≤ δ ≤ 20 (q = 75). Proposition 3.3(i). On the other hand, for fixed (q, α) = (75, 1), the stability and instability regions with varied δ are shown in Figure 6(b). We can see from this figure that E * is stabilized as δ becomes large, which is consistent with Proposition 3.3(ii).
Finally, to confirm the validity of Proposition 4.1, we fix (q, α) = (75, 1) and vary L to illustrate the stability and instability regions determined by T + n and T − n , n = 0, 1, 2 (see Figure 7).
We can confirm from Figure 7 that E * is asymptotically stable when T and L are small. This result is consistent with the statement in Proposition 4.1.

Discussion
In this paper, we have constructed an epidemic model with quarantine and distributed time delay. We have defined the basic reproduction number R 0 and shown that if R 0 ≤ 1, then the disease-free equilibrium E 0 is globally asymptotically stable, whereas if R 0 > 1, then E 0 is unstable and a unique endemic equilibrium E * exists. This implies that the backward bifurcation at R 0 = 1 never occurs in our model and reducing R 0 less than 1 is a clear goal to control the epidemic. Moreover, we have investigated the local asymptotic stability of E * for R 0 > 1, and shown that either if the sensitivity q of the quarantine is sufficiently small or the lifting rate δ is sufficiently large, then E * is locally asymptotically stable. These conditions correspond to the situation where the quarantine is less restrictive. On the other hand, we have obtained sufficient conditions for a Hopf bifurcation and shown by numerical simulations that periodic solutions can exist in some parameter regions. From Figure  6, one can see that unstable parameter regions are enlarged as q increases or δ decreases. This observation suggests that recurrent epidemic waves become likely to happen when the quarantine is meaningful and the time delay effect exists. In fact, Propositions 3.3 and 4.1 ensure that periodic solutions do not bifurcate from the endemic equilibrium when there are less quarantine and time delay effects. In conclusion, our results suggest that the quarantine and time delay could play important roles in the occurrence of recurrent epidemic waves.
In this study, we have chosen the time delay T as the bifurcation parameter. For practical purposes, bifurcation analysis with respect to control parameters q and δ may be more significant. Actually, Figure 6 reveals that the stability change can occur for fixed T by changing q and δ (for example, if we fix T = 5 in Figure 6(a), the stability change occurs at q ≈ 75). In other words, when the time delay T is fixed for a given disease, we can estimate the critical values of control parameters q and δ that induce recurrent epidemic waves.
Note that parameters T, q and δ do not affect on the basic reproduction number R 0 . That is, the time delay and quarantine in our model do not affect on the average epidemic size, but can determine the existence of recurrent epidemic waves. Our results suggest that if we want to avoid recurrent epidemic waves, it could be better to make the quarantine less restrictive and let the disease attain a steady state early. This strategy could contribute to the early settlement of the problem of the coexistence of health and economy, which has been a central problem raised in COVID-19. However, of course, restrictive quarantine measures would be effective in the initial invasion phase of a new epidemic to gain time for improving the medical system and developing vaccines.
Finally, we state some limitations and open questions in our work. In this paper, to make the analysis simple, we have disregarded some important structures such as seasonality, viral mutation, age and space. Incorporating such structures into our model would be an important future work. In particular, we have restricted our attention to the spread of a disease on a homogeneous network. Considering the heterogeneous network would lead to more complicated dynamics. Such consideration would be needed if we perform a more detailed analysis on the sensitivity of the heterogeneity to the epidemic dynamics. In addition, we have not used any real data in the numerical simulations. Applying our model to a specific disease such as COVID-19 would also be an important future work. A such attempt could give a new insight on the validity of our conclusion in this paper.