A MULTI-GROUP SIR EPIDEMIC MODEL WITH AGE STRUCTURE

. This paper provides the ﬁrst detailed analysis of a multi-group SIR epidemic model with age structure, which is given by a nonlinear system of 3 n partial diﬀerential equations. The basic reproduction number R 0 is obtained as the spectral radius of the next generation operator, and it is shown that if R 0 < 1, then the disease-free equilibrium is globally asymptotically stable, while if R 0 > 1, then an endemic equilibrium exists. The global asymptotic stability of the endemic equilibrium is also shown under additional assumptions such that the transmission coeﬃcient is independent from the age of infective individuals and the mortality and removal rates are constant. To our knowledge, this is the ﬁrst paper which applies the method of Lyapunov functional and graph theory to a multi-dimensional PDE system.


1.
Introduction. It is well-known that the spread of infectious diseases can be modeled by differential equations. The SIR model is one of the most popular epidemic models, in which total host population is divided into three classes called susceptible (S), infective (I) and removed (R). For general information of mathematical epidemiology, the reader may refer to Diekmann et al. [6].
The age-structured SIR models have been studied by many authors (see [1,4,7,10,11,15,16,18,28,31,35,40,41,42,46]). A typical age-structured SIR model is given by the following nonlinear system of partial differential equations: where S (t, a), I (t, a) and R (t, a) denote the densities of susceptible, infective and removed individuals of age a at time t, respectively. µ (a) denotes the age-specific per capita mortality rate, γ (a) denotes the age-specific per capita removal rate and f (a) denotes the age-specific per capita fertility rate. The density of total host population of age a at time t is given by P (t, a) := S (t, a) + I (t, a) + R (t, a). The force of infection to a susceptible individual of age a at time t is given by where β (a, σ) denotes the disease transmission coefficient between a susceptible individual of age a and an infective individual of age σ. For model (1.1), the basic reproduction number R 0 is obtained as the spectral radius of a linear operator called the next generation operator (see Diekmann et al. [6]). The global asymptotic stability of the disease-free equilibrium of (1.1) was proved for R 0 < 1 in Inaba [18]. However, the global asymptotic stability of an endemic equilibrium of (1.1) for R 0 > 1 has been an open problem for a long time and even the possibility of unstable endemic equilibria has been shown for R 0 > 1 (see Thieme [40], Andreasen [1], Cha et al. [4], Franceschetti et al. [9]). On the other hand, in Melnik and Korobeinikov [35] the global asymptotic stability of the endemic equilibrium for systems in which only S has the age structure was proved by constructing a Lyapunov functional. Under the assumption that host population is not homogeneous but heterogeneous, we can divide each of the classes (S, I and R) into several homogeneous groups according to the heterogeneity (e.g., sex, position) of each individual. Such heterogeneity in the host population can result from different contact modes such as those among children and adults for childhood diseases (e.g., measles and mumps), or different behaviors such as the numbers of sexual partners for some sexually transmitted infections (e.g., Herpes, Condyloma acuminatum). Such models are typically called multi-group epidemic models (see van den Driessche and Watmough [43]). The multi-group SIR epidemic models as ODEs have been studied by many authors (see [2,12,14,45,25,28,39]). For a multi-group SIR epidemic model, Guo et al. [12] proved the complete global asymptotic stability of each equilibrium by using a graph-theoretic approach. By using the results or ideas of Guo et al. [12], the papers [39,45,28] investigated uniqueness and global stability of the endemic equilibrium for several classes of multigroup models, and some open problems were resolved.
In this paper, we are concerned with an SIR epidemic model in quite a general form, containing both of the age structure and the multi-group structure. This model is suitable for a disease which exhibits the age dependency and the heterogeneous state dependency. An SIR epidemic model we shall consider in this paper is formulated by the following nonlinear system of first-order partial differential equations: (1.2) where S j (t, a), I j (t, a) and R j (t, a) denote the densities of susceptible, infective and removed individuals of age a at time t in group j, respectively. µ j (a) and γ j (a) denote the age-specific per capita mortality and removal rates for individuals in group j, respectively. f jk (a) denotes the age-specific per capita fertility rate at which an individual of age a in group k gives birth to an offspring in group j. P j (t, a) := S j (t, a) + I j (t, a) + R j (t, a) denotes the density of total host population of age a at time t in group j, and λ j (t, a) = n k=1 +∞ 0 β jk (a, σ) I k (t, σ) dσ, j = 1, 2, . . . , n denotes the force of infection to a susceptible individual of age a at time t in group j, where β jk (a, σ) denotes the disease transmission coefficient between a susceptible individual of age a in group j and an infective individual of age σ in group k.
To our knowledge, there are very few studies on epidemic models with both the age-structure and multi-group structure. In Feng et al. [8], the global behavior of an SIS epidemic model with both of these structures was studied. In Li et al. [30], an age-structured epidemic with two-group structure including super-infection was studied. An SIR epidemic model similar to (1.2) was studied in Kuniya [28], where the main focus was put on a corresponding ODE model derived by using a discretization approach as in Tudor [41] and Hethcote [15]. Therefore, a detailed mathematical analysis for the PDE model (1.2) has been an open problem. By virtue of both of the age-structure and multi-group structure, model (1.2) is thought to be realistic for sexually transmitted diseases, in which the age of sexual activity and the heterogeneity as sex play important roles. In this viewpoint, the detailed mathematical analysis for model (1.2) is an important task not only for mathematical but also biological purposes. In this paper, we shall investigate the existence, uniqueness and asymptotic stability of each equilibrium of (1.2) by using a functional analytic approach as in Webb [46] and Inaba [18] and a Lyapunov functional approach as in Magal et al. [31], McCluskey [33] and Melnik and Korobeinikov [35].
The organization of this paper is as follows. In Section 2, we prepare some preliminaries. In Section 3, using an abstract formulation, we show the mathematical well-posedness of the problem. In Section 4, we derive the basic reproduction number R 0 for model (1.2) as the spectral radius of the next generation operator. In Section 5, we investigate the existence of each equilibrium. It is shown that the basic reproduction number R 0 plays the role of a threshold value for the existence of a nontrivial endemic equilibrium. In Section 6, we investigate the uniqueness of the endemic equilibrium for two cases of separable mixing and general mixing. In Section 7, we prove the global asymptotic stability of the disease-free equilibrium for R 0 < 1. In Section 8, under the assumptions that the disease transmission coefficients β jk (a, σ) are independent of the age of infective individuals and the mortality and removed rates are constant, we prove the global asymptotic stability of an endemic equilibrium for R 0 > 1. In Section 9, we perform numerical simulation to verify the validity of our main theorem.

Preliminaries.
We make the following assumption on each coefficient of system (1.2): From (i) and (ii) of Assumption 2.1, we see that for each j, k ∈ {1, 2, . . . , n}, there exist positive constants γ + j > 0 and β + jk > 0 such that ess.sup By adding the first three equations of (1.2), we obtain the following Lotka-McKendrick-von Foerster system: Note that from (ii) of Assumption 2.1, the interval of integration in the second equation of (2.2) is bounded above by a † . Let and M (a) := (f jk (a)) 1≤j,k≤n . We make the following assumption: 2) has a nonnegative demographic steady state P * j ∈ L 1 + (0, +∞)\{0} , j = 1, 2, . . . , n, which is globally stable. In what follows, we assume that the demographic steady state has been already reached, that is, P j (t, ·) ≡ P * j (·) , j = 1, 2, . . . , n. Thus, system (1.2) can be rewritten as follows: are positive constants.
3. Well-posedness of the problem. To define the state space for (2.4), we prepare the function spaces respectively, where T denotes the transpose of a vector. Let us denote their positive cones by X + , Y + and Z + , respectively. Under the assumption that the system has reached the demographic steady state, we define the following state space: where S := (S 1 , . . . , S n ) T ∈ X, I := (I 1 , . . . , I n ) T ∈ X, R := (R 1 , . . . , R n ) T ∈ X and P * := (P * 1 , . . . , P * n ) T ∈ X. In what follows, we assume that the initial condition of (2.4) is taken from Ω.
To show the existence and uniqueness of a global classical solution of problem (2.4) in Ω, it is sufficient to consider the following reduced 2n-dimensional system: Note that S j (j = 1, 2, . . . , n) is given by S j = P * j − I j − R j in Ω. Let us define the state space for (3.1) bỹ Let us define a linear operator A : ∈ Y : ψ 1 and ψ 2 are absolutely continuous,  Then, setting u := (I, R) T ∈ Y , we rewrite problem (3.1) as the following abstract Cauchy problem inΩ: where u 0 = (I 1,0 , . . . , I n,0 , R 1,0 , . . . , R n,0 ) T ∈ Y . We easily see that operator A is an infinitesimal generator of a C 0 -semigroup e tA t≥0 : Y → Y defined by From (3.8), the spaceΩ is positively invariant under the semiflow defined by e tA , that is, e tA Ω ⊂Ω for all t ≥ 0. where I denotes the identity operator.
Let α > 0 be a constant such that (3.9) holds. Using this α, as in Busenberg et al. [3], we rewrite the abstract Cauchy problem (3.7) as and investigate the existence of solution of problem (3.11). The mild solution of (3.11) is given by the solution of integral equation As in Busenberg et al. [3], we construct an iterative sequence {u n : n ∈ N ∪ {0}} such that From Lemma 3.1, and by using the facts that u n+1 is a linear convex combination of e tA u 0 ∈Ω and e (t−s)A (I + αF ) u n ∈Ω, we conclude that u n+1 ∈Ω if u n ∈Ω. Thus, from the Lipschitz continuity of the operator F , u n converges to the mild solution u ∈Ω of (3.11) uniformly as n → +∞. It is easy to see that if the mild solution u ∈Ω of (3.11) belongs to D (A), then it is differentiable and hence it is a classical solution of problem (3.7). In conclusion, we have the following proposition: 4. The basic reproduction number R 0 . In this section, we define the basic reproduction number R 0 , which is epidemiologically the expected number of secondary cases produced by a typical infective individual during its entire period of infectiousness in a fully susceptible population. The basic reproduction number is the most important idea characterizing the renewal process in structured populations, so it has been developing as a central dogma in mathematical epidemiology since an epoch-making paper by Diekmann, Heesterbeek and Metz at 1990 ( [5]). For recent developments of R 0 , the reader may refer to Diekmann et al. [6] and Inaba [22].
where Q (a) and Γ (a) are diagonal matrices defined by (2.3) and (3.6), respectively. Then, the linear system (4.1) is rewritten as the following abstract Cauchy problem in X: where I = (I 1 , . . . , I n ) T ∈ X, I 0 = (I 1,0 , . . . , I n,0 ) T ∈ X and P is a bounded linear operator on X defined by Let V (t) = e tB be the C 0 -semigroup generated by the generator B. Using the variation of constants formula, we have Applying P from the left hand side to (4.5), we obtain an abstract renewal equation: where v(t) := P I(t) denotes the density of newly infecteds in the linear invasion phase, g(t) := P V (t)I 0 and Ψ(s) := P V (s). Then the next generation operator is defined by where we have used a relation that for z ∈ ρ(B) (ρ(B) denotes the resolvent set of B) and 0 ∈ ρ(B).
Then the next generation operator is calculated as follows: (4.9) According to Diekmann, Heesterbeek and Metz [5], we define the basic reproduction number R 0 by the spectral radius of the next generation operator K: where r(A) denotes the spectral radius of a bounded linear operator A.
In fact, we can show that R 0 gives the asymptotic per generation growth factor of infected population in the linear invasion phase, the Malthusian parameter (asymptotic growth rate of infecteds) is positive R 0 > 1, while it is negative if R 0 < 1 (see Inaba [22]).

5.
Existence of equilibria. Instead of system (2.4), we next focus on the following reduced 2n-dimensional system: In this section, we investigate the existence of equilibria of system (5.1).
Denote S * := (S * 1 , . . . , S * n ) T ∈ X and I * := (I * 1 , . . . , I * n ) T ∈ X the equilibria of system (5.1). They must satisfy the following equations: It is obvious that S * j (a) ≡ P * j (a) , I * j (a) ≡ 0, j = 1, 2, . . . , n is a solution set of (5.2) and hence, it is the disease-free equilibrium of system (5.1). From the variation of constants formula, we have and Substituting (5.4) into the third equation of (5.2), we have (5.6) Set the following linear operator on X corresponding to the right-hand side of this equation: where Then, λ * j = H j (λ * ), j = 1, 2, · · · , n. From (5.3) and (5.5), we havẽ Hence, defining a nonlinear operator we can show the existence of the endemic equilibrium by finding a positive nontrivial fixed-point of operator Φ since it implies the existence of positiveλ * . Noticing that H is linear and H(0) = 0, we see that the Fréchet derivative of Φ at ϕ = 0 is Noticing that we see that Φ [0] is equal to the next generation operator K defined by (4.7). In what follows, we shall show that the basic reproduction number R 0 = r(K) determines the existence of the positive fixed point of operator Φ. Before the proof of our main result in this section, we adopt the following assumption: . . , n}, the following holds uniformly for σ ≥ 0: Under Assumption 5.1, we first prove the following lemma: Proof. Let B 0 be an arbitrary bounded subset of X. Note that there exists a positive n is continuous, monotone nonincreasing and lim a→+∞ P * j (a) = 0, Assumption 5.1 implies that for any > 0, there exists a δ > 0 such that holds whenever |h| < δ. Hence, from (5.10) and Assumption 2.1 (i), the inequality (5.9) can be further evaluated as Since the right-hand side of this inequality converses to zero as A → +∞ uniformly for ϕ ∈ B 0 , the assertion is true. Consequently, from the well-known compactness criterion for L 1 -space (see, e.g., Yosida [47, p.275]), we conclude that operator K is compact.
Next we prove (ii). Let X * be the dual space of X, that is, the set of all linear functionals on X. Let f, ϕ denote the value of f ∈ X * at ϕ ∈ X. Let X * + be the dual cone of X, that is, the subset of X * consisting of all positive linear functionals on X. We write f ∈ X * + if f ∈ X * and f, ϕ ≥ 0 for all ϕ ∈ X + .
Next we prove the following lemma.
Then, Φ is given by the compositionΦ • H. To show (ii), it suffices to show that Φ(X + ) is bounded since H(X + ) ⊂ X + . In fact, which implies thatΦ(X + ) is bounded and hence, Φ(X + ) is bounded.
To show (i), it suffices to show that H is compact since the composition of a bounded operator and a compact operator is also compact. In fact, under Assumption 5.1 (ii), we can show the compactness of H in a similar way as in the proof of Lemma 5.1 (i).
In order to show the existence of a nontrivial fixed-point of Φ, we shall apply the Krasnoselskii fixed-point theorem (see Krasnoselskii [26,Theorem 4.11]). Before the proof, we recall some basic facts on the theory of linear positive operators. Let E be a real or complex Banach space.
The cone C is called solid if its interior C 0 is nonempty. Combining the Krein-Rutman theorem (see Krein and Rutman [27]) and Sawashima's theorem (see Sawashima [37]), we can obtain the following corollary (see Inaba [21,Corollary 7.6]): Corollary 5.1. Let C be a cone of a Banach space E and T be a bounded linear operator from E to itself. Suppose that C is total, the spectral radius r (T ) > 0 of operator T is positive and T is compact and nonsupporting with respect to C. Then, the following holds: (i) r (T ) ∈ P σ (T ) \ {0}, where P σ (·) denotes the point spectrum of an operator. Moreover, r (T ) is a simple pole of the resolvent R (λ, The eigenspace corresponding to r (T ) is one-dimensional and its eigenvector ψ ∈ C is a quasi-interior point, that is, f, ψ > 0 for all f ∈ C * \ {0}.
Moreover, any eigenvector of T in C is proportional to ψ. (iii) The adjoint eigenspace corresponding to r (T ) is one-dimensional and its eigen- −n T n converges in the topology of the operator norm and it is a strictly nonsupporting operator given by where Γ 0 denotes the positively oriented circle centered at r (T ) such that no points of the spectrum σ (T ) except r (T ) lie on and inside Γ 0 .
As a consequence of the Krasnoselskii fixed-point theorem, we can obtain the following corollary (see also Inaba  Proof. For nonlinear operator Φ and real Banach space X, Corollary 5.2 (i) follows from Lemma 5.2. Corollary 5.2 (ii) immediately follows with the strong Fréchet derivative K = Φ [0]. Since the compactness and nonsupporting property of K follows from Lemma 5.1, Corollary 5.1 can be applied to show that r (K) (= R 0 > 1) is the only positive eigenvalue of K corresponding to a positive eigenvector and there is no eigenvector of K corresponding to eigenvalue 1. Therefore, Corollary 5.2 (iii) follows. Consequently, Φ has at least one nontrivial fixed-point in X + \ {0}.

6.2.
Case of general mixing. Next we consider the case of general mixing. In this section, we make the following assumption: Assumption 6.2. For each j, k ∈ {1, 2, . . . , n},
Under Assumption 6.2, we first prove the following lemma: Lemma 6.1. LetΦ be defined by (6.3). Under Assumption 6.2,Φ is monotone nondecreasing on X.

7.
Global stability of the disease-free equilibrium. In this section, we investigate the global asymptotic stability of the disease-free equilibrium (P * , 0) ∈ Y + of system (5.1).
Using the generator B introduced in Section 4, the second equation of system (5.1) can be rewritten as the following abstract Cauchy problem in X:    d dt I (t) = BI (t) + (ΛI (t)) S (t) , where S = (S 1 , . . . , S n ) T ∈ X, I = (I 1 , . . . , I n ) T ∈ X, I 0 = (I 1,0 , . . . , I n,0 ) T ∈ X and Λ is a linear operator on X defined by (3.5). If we regard S (t) as a given function, we can regard (7.1) as a linear non-autonomous equation for I (t). Since S (t) ≤ P * , we have From Proposition 4.13 in Webb [46], we have where σ (A) and E σ (A) denote the spectrum and the essential spectrum of a closed linear operator A, respectively. To prove the global asymptotic stability of the disease-free equilibrium, it is sufficient to show that the growth bound ω 0 (B + P ) is negative. To this end, we first prove the following lemma: Proof. Observe that , a − t > 0. As in the proof of Theorem 4.6 in Webb [46], we have (7.6) Moreover, we can prove the compactness of operator P based on the assumption 5.1 and the fact that P * j ∈ L 1 ∩ L ∞ , so it follows that ω 1 (B + P ) = ω 1 (B) (see Webb [46,Proposition 4.14]). Hence we obtain (7.5).
Note that U 0 = K, where K is the next generation operator defined by (4.8).
Proof. First note that the compactness of operator U z can be proved as in the proof of (i) of Lemma 5.1. Next let us prove that if z ∈ Σ, then z ∈ P σ (B + P ) ∩ z ∈ C : z > −µ . Suppose that there exists an eigenvector which implies that Pφ is an eigenfunction associated with eigenvalue unity of U z , that is, z ∈ Σ.
Therefore ω 0 (B + P ) > 0 if R 0 > 1 and ω 0 (B + P ) = 0 if R 0 = 1. On the other hand, if R 0 < 1, there is no normal eigenvalue with nonnegative real part. Since there exist at most finite number of normal eigenvalue in the half plane z > −µ, it follows from (7.4) that ω 0 (B + P ) < 0.
Next we prove the following lemma.   We show (iii). Since we shall consider h → +0, we can assume that h ∈ (0, t) without loss of generality. Then, from (8.12), we have where β + j (j = 1, 2, . . . , n) is a positive upper bound of β j (a) whose existence is assured by Assumption 2.1. Note that it follows from the second equation of (8.3) and therefore, J (t) is Lipschitz continuous on R + with Lipschitz coefficient M J . Hence, from (8.14), we obtain Consequently, it follows from (8.13), (8.15) and (8.16) that → 0 as h → +0 uniformly in x 0 ∈ C, and hence, (iii) is shown and the proof is complete.
Using Lemmas 8.1 and 8.2, we prove the following proposition on the relative compactness of the positive orbit {U (t) x 0 : t ≥ 0}: Proposition 8.1. Under Assumption 8.1, let U (t) and C be defined by (8.5) and (8.6), respectively. For x 0 ∈ C, {U (t) x 0 : t ≥ 0} has compact closure in Y.
for any x 0 ∈ C, it follows from Lemmas 8.1-8.2 and Proposition 3.13 in Webb [46] that {U (t) x 0 : t ≥ 0} has compact closure in Y.
Since the relative compactness of the solution orbit is proved in Proposition 8.1, we can use the invariance principle stated in Walker [44,Theorem 4.2 in Chapter IV] together with a Lyapunov functional in order to prove the global asymptotic stability of endemic equilibrium (S * , J * ). Thus, following the graph-theoretic approach of Guo et al. [12], we construct a Lyapunov functional. Set constants β jk := +∞ 0 β j (a) S * j (a) da J * k , j, k = 1, 2, . . . , n and a Laplacian matrix Then, from Guo et al. Using this κ, we consider the following Lyapunov functional: where g (w) := w − 1 − ln w. To check that this function L is well-defined, it suffices to show that the integration in the right-hand side is bounded. We make the following assumption: For each j ∈ {1, 2, · · · , n}, S j,0 ∈ X + , J j (0) > 0. In addition, For a − t ≥ 0, b k µ k + µ j ae −µj a , j = 1, 2, · · · , n.
Under Assumption 8.2, the second term in the right-hand side is finite. Thus, S * j (a) ln S j (t.a)/S * j (a) , j = 1, 2, · · · , n is integrable and hence, the Lyapunov functional L(S, J) is well-defined.
Under these preparations, we prove the following proposition: Proof. The uniqueness of (S * , J * ) follows from Proposition 6.1. To prove the global stability of (S * , J * ), we use the Lyapunov functional L (S, J) defined by (8.19). The differentiation of L along the trajectories of system (8.3) is  (iv) If R 0 > 1 and Assumption 8.1 holds, then system (5.1) has a unique globally asymptotically stable endemic equilibrium (S * , I * ) ∈ Y + \ {0} which attracts any initial values satisfying Assumption 8.2.

9.
Numerical examples. In this section, we verify the validity of Theorem 8.1 by performing numerical simulation. Since we can not consider the infinitely large age interval a ∈ [0, +∞) in numerical simulation, we set the maximum age a M and normalize it as a M = 1 for simplicity. Let us consider two-group case (n = 2) for sexually transmitted diseases, in which j = 1 and j = 2 imply female and male populations, respectively. Let us fix the following demographic parameters: b 1 = 1.5, b 2 = 1, µ 1 (a) ≡ µ 1 = 5, µ 2 (a) ≡ µ 2 = 10. (9.1) In this case, the demographic steady states P * j (a) = b j exp (−µ j a), j = 1, 2 are almost zero at a = a M = 1 (see Figure 1). Therefore, we can expect that our simulation would be a good approximation for the case of infinite age interval [0, +∞) in which P * j (a) → 0, j = 1, 2 as a → +∞. Let us fix the following epidemic parameters γ 1 = 0.25, γ 2 = 0.5, β 21 (a, σ) = β 22 (a, σ) = β 2 (a) = 5 1 + 0.9 cos 3πa 2 and vary β 11 (a, σ) = β 12 (a, σ) = β 1 (a). Note that Assumption 8.1 holds in this case. Therefore, from Theorem 8.1, we can expect that the basic reproduction number R 0 would play the role of the complete threshold value for the global asymptotic stability of each equilibrium. From (4.9), the component K j ϕ in the vector of the next generation operator (4.8) can be approximated by the following one.
On the other hand, for β 1 (a) = 11.5 (1 + 0.9 cos 3πa/2), we have R 0 ≈ 1.0465 > 1. Hence, from Theorem 8.1 (iv), we see that the endemic equilibrium is globally asymptotically stable. In fact, in Figure 3, both of the infective individuals converge to the positive distributions.
10. Discussion. In this paper, we have formulated a multi-group SIR epidemic model (1.2) with age structure. For the reduced system (5.1), we have shown that the basic reproduction number R 0 defined by (4.10) plays the role of a threshold value for the existence, uniqueness and asymptotic stability of each equilibrium. . Time variation of I 1 (t, a) and I 2 (t, a) for β 1 (a) = 11.5 (1 + 0.9 cos 3πa/2) (in this case, R 0 ≈ 1.0465 > 1) In this paper, we have mainly focused on the stability of each equilibrium. However, in other references, some irregular cases such as unstable endemic equilibrium (see Andreasen [1], Cha et al. [4] and Thieme [40]) or multiple endemic equilibria (see Franceschetti et al. [9]) when R 0 > 1 have been observed for age-structured SIR epidemic models. Investigating such nontrivial properties of our model (1.2) has been left as an open problem.