Modelling diseases with relapse and nonlinear incidence of infection: a multi-group epidemic model

In this paper, we introduce a basic reproduction number for a multi-group SIR model with general relapse distribution and nonlinear incidence rate. We find that basic reproduction number plays the role of a key threshold in establishing the global dynamics of the model. By means of appropriate Lyapunov functionals, a subtle grouping technique in estimating the derivatives of Lyapunov functionals guided by graph-theoretical approach and LaSalle invariance principle, it is proven that if it is less than or equal to one, the disease-free equilibrium is globally stable and the disease dies out; whereas if it is larger than one, some sufficient condition is obtained in ensuring that there is a unique endemic equilibrium which is globally stable and thus the disease persists in the population. Furthermore, our results suggest that general relapse distribution are not the reason of sustained oscillations. Biologically, our model might be realistic for sexually transmitted diseases, such as Herpes, Condyloma acuminatum, etc.


Introduction
In recent years, great attention has been paid in analysing the multi-group epidemic models, which have been proposed to describe the disease transmission dynamics of many infectious diseases in heterogeneous host populations, such as measles, mumps and gonorrhoea and vectorborne diseases such as West-Nile virus and Malaria. For more and detailed justifications for multi-group disease models and many different types of heterogeneity epidemic models (see e.g. [3,13,14,16,20,23,32] and references cited therein). It is well known that long-time behaviours of multi-group models with higher dimensions, especially the global asymptotical stability of the endemic equilibrium (EE), is a very challenging topic. The question of uniqueness and global stability of the EE, when the basic reproduction number is more than one, has largely been open. In [13], a graph-theoretic approach was developed, which is known to be an effective tool for the global stability analysis of multi-group epidemic models. Much research has been done with this approach [3,13,14,22,23,35], focusing on understanding the transmission mechanism and the global behaviour of multi-group epidemic models. So the study of relapse distribution and different nonlinear incidence of infection have been the subject of intense theoretical analysis on the heterogeneous epidemic models in the literature [26,34,35].
In this paper, motivated by the works of van den Driessche and Zou [10], Wang et al. [34] and Sun and Shi [31], we shall investigate the global dynamics of a general multi-group epidemic model with general relapse distribution and nonlinear incidence rate, in particular to investigate the impacts of heterogeneity and nonlinear incidence rate on the dynamics of the basic SIR epidemic model. The population is divided into n distinct groups (n ≥ 1). For 1 ≤ i ≤ n, the ith group is further partitioned into three compartments: The class of ϕ i (S i ) that satisfy (A 1 ) include both λ i − d S i S i and λ i − d S i S i + r i S i (1 − S i /N i ), which have been widely used in the literature of population dynamics [1,13].
Since nonlinear incidence of infection has been observed in disease transmission dynamics, it has been suggested that the standard bilinear incidence rate shall be modified into a nonlinear incidence rate in many research [18,28,29]. In this paper, we replace the incidence rate f (S)I in [34] by a general form f (S, I). We assume that the disease incidence in the ith group can be calculated as where the sum takes into account cross-infections from all groups, and β ij represents the transmission coefficient between compartments S i and I j . Throughout the paper, β ij is non-negative for all i, j = 1, 2, . . . , n, and n-square matrix B = (β ij ) 1≤i,j≤n is irreducible [4]. Biologically, this is the same as assuming that any two groups i and j have a direct or indirect route of transmission. More specifically, individuals in I j can infect ones in S i directly or indirectly.
Denote P i (t) by the fraction of recovered individuals remaining in the recovered class t time units after recovery in each group. It was assumed in [10] that P i (t) satisfies the following reasonable properties: is non-increasing, piecewise continuous with possibly finitely many jumps and satisfies P i (0 dt is positive and finite.

101
The proportion of recovered individuals in the ith group can be expressed by the integral where γ i > 0 is the recovery rate constant assuming that the infective period is exponentially distributed in the ith group. The term e −d i (t−ξ) in the above integral accounts for the death of infectives in ith group. It is assumed that no individuals are initially in the recovered class, i.e.
Here, the integral is in the Riemann-Stieltjes sense to allow for possible jump discontinuities of P i (t).
Hence, the new multi-group epidemic model with relapse distribution and nonlinear incidence rates can be written as the following 3n-dimensional system of differential and integral equations: which is an improvement of our earlier model studied in [34]. It also covers many research works in the literature, for example, the ones in [10,26]. The disease transmission diagram is depicted in Figure 1. And the parameters in the model are summarized in the following list: β ij : coefficient of transmission between compartments S i and I j ; d i : natural death rates of all compartments in the ith group; γ i : rate of recovery of infectious individuals in the ith group.
In the present paper, our goal is to carry out a complete mathematical analysis of system (1) and establish its global dynamics. By virtue of both the relapse distribution and multi-group structure, model (1) is thought to be realistic for sexually transmitted diseases such as herpes simplex virus type 2 (Herpes) [6,7,33], for which disease is transmitted by close physical or sexual contact, recovered individuals may revert back to the infective class because of reactivation of the latent infection or incomplete treatment [7,33]. One important feature of herpes is that an individual once infected remains infected for life, and the virus reactivates regularly with reactivation producing a relapse period of infectiousness (see e.g. [6] and references cited therein). Therefore, the detailed mathematical analysis for model (1) is an important task from not only mathematical but also biological points. It is thus of interest to investigate whether sustained oscillations are the result of general relapse distribution. This provides us with one motivation to conduct our work.
The paper is organized as follows. In the next section, we give some preliminaries on system (1). Our main results are stated in Section 3. One case for n = 1 is investigated in Section 4. Finally, we give a brief summary and discussion.

Preliminaries
We first consider system (1) in the phase space R 3n and with the initial conditions On the basis of biological considerations, we make the following assumptions for f ij (S i , I j ) [31]: Typical examples of f ij (S i , I j ) satisfying (A 2 )-(A 4 ) include common incidence functions such as f ij (S i , I j ) = S i I j [12,13,19]; f ij (S i , I j ) = ηS i I j /(1 + θ S i ) [1,11]; f ij (S i , I j ) = S q i I j [24,25,37]. (2) with S i (0) ∈ R + , I i (0) ∈ R + and R i (0) ∈ R + , the solutions of system (1) are ultimately uniformly bounded in R 3n .

Lemma 2.1 For initial conditions in Equation
Proof It follows from (A 1 ) and the first equation of Equation (5) that lim sup t→∞ S i (t) ≤ S 0 i for all i = 1, 2, . . . , n. Next we show that the solution of system (1) is ultimately bounded in R 3n + . For each i, adding the three equations in Equation (1) gives where (3), we can also obtain that, for any is a forward invariant compact absorbing set with respect to the system (1). Let It can be shown that 0 is the interior of . Furthermore, all positive semi-orbits in are precompact R 3n [2] and thus have non-empty ω-limit sets.
Note that the first two equations of system (1) are independent from R i , and therefore, the dynamics is governed by the reduced system The initial condition of system (5) is assumed to be given as For system (5) with properties (H 1 ), the existence, uniqueness and continuity of solutions follow from the theory for integro-differential equations in [30]. Moreover, it can be verified that every solution of Equation (5) with non-negative initial data remains non-negative. It follows from Equation (4) that the set is the subset of . Thus, is also the subset of 0 . Furthermore, it is a forward invariant compact absorbing set with respect to the system (5). All positive semi-orbits in 0 are precompact R 2n and thus have non-empty ω-limit sets. System (5) always admits a disease-free equilibrium (DFE), P 0 = (S 0 1 , 0, . . . , S 0 n , 0), in 0 . Let Clearly,P i are the average time that recovered individuals remain in the recovered class before relapsing. By the properties of P i , one knows that Actually, d iPi are the probability that recovered individuals will die during the recovery period. Hence, Q i represent the proportion of the recovered individuals who could survive the recovery period, where The basic reproduction number R 0 is defined as the expected number of secondary cases produced in an entirely susceptible population by a typical infected individual during their entire infectious period [8,9]. For system (5), we can calculate it as the spectral radius of a matrix called the next generation matrix. Let then the next generation matrix is and hence the basic reproduction number of model (5) is calculated by the spectral radius of the next generation matrix where ρ(·) denotes the spectral radius of matrix. It is well known that ρ( where Since it can be verified that system (5) satisfies conditions (A 1 )-(A 5 ) of [9, Theorem 2], we have the following proposition.
Note that Equation (5) may not have an EE for finite time t. It follows from [30] that if Equation (5) has an EE, then the EE must satisfy the limiting system given by Since the limiting system (9) contains an infinite delay, its associated initial condition needs to be restricted in an appropriate fading memory space. For any λ i ∈ (0, d i ), define the following Banach space of fading memory type (see e.g. [15,17] and references therein) Standard theory of functional differential equations [17] implies I it ∈ C i for t > 0. Thus, we will consider system (9) in the phase space It can be verified that solutions of Equation (9) in with initial conditions (11) remain nonnegative. An equilibrium P * = (S * 1 , I * 1 , . . . , S * n , I * n ) in the interior of 0 is called an EE, where S * i , I * i > 0 satisfy the equilibrium equations

Global stability results
The global dynamical behaviour of system (5) and (9) is completely established in the following results. Denote Obviously, H : R + → R + attains its strict global minimum at z = 1 and H(1) = 0.

Global dynamics of DFE
We first claim that there does not exist any EE P * in . If we assume that S = S 0 , it follows Let Consider a Lyapunov functional where U + is given as It is easy to know that L DFE ≥ 0 with equality if and only if I i (t) = 0 and J i (ξ )I i (t − ξ) = 0 for almost all ξ ≥ 0.
Next we show that the EE P * of system (9) is unique and globally asymptotically stable when R 0 > 1 . Summarizing the statements, uniform persistence of Equation (9) from Theorem 3.1, together with uniform boundedness of solutions in the interior of , implies system (9) admits at least one EE [5, Theorem 2.8.6].
Let X(t) = (S i (t), φ it ) be a solution of Equation (9). By Theorem 3.1 and using similar arguments to [27], it follows that the ω-limit set W of X is non-empty, compact and invariant and that W is the union of orbits of Equation (9).

Global dynamics of EE
To get the global stability of P * , we make the following assumptions [31]: For convenience of notations, set Then,B is also irreducible. One knows that the solution space of the linear system Proof Let P * = (S * 1 , I * 1 , . . . , S * n , I * n ) denote the unique EE of system (9). Define a Lyapunov functional as The definition of the fading memory space implies that L EE is well-defined, that is, L EE is bounded for all t ≥ 0. It follows from Lemma 2.1 of [34] and assumptions ( Differentiating L S along the solution of system (9) and using equilibrium equations (12) and (13), we obtain dL S dt (9) = Differentiating L I along the solution of system (9), we obtain Differentiating L − along the solution of system (9) and using integration by parts, we obtain Combining Equations (16)- (18) yields Then, by assumption (A 5 ) − (A 6 ), Furthermore, under (A 7 ), we have where G i (I i ) = −I i /I * i + ln(I i /I * i ) Obviously, the equalities in Equations (19) and (20) hold if and only if i.e. S i = S * i , I i = I * i , i = 1, 2, . . . , n. We can show that F ij andβ ij satisfy the assumptions of Theorem 3.1 and Corollary 3.3 in [22]. Therefore, the function L = n i=1 v i (L S + L I + L − ) is a Lyapunov function for system (9), namely, L EE | (9) ≤ 0 for all (S 1 , I 1 , S 2 , I 2 , . . . , S n , I n ) ∈ . One can only show that the largest invariant subset where L EE | (9) = 0 is the singleton P * using the same argument as in [22]. By LaSalle's invariance principle, P * is globally asymptotically stable in 0 . This completes the proof of Theorem 3.2.

Special case for n = 1
When n = 1, system (5) reduces to a single-group SIR model with general relapse distribution and nonlinear incidence rate: Following the method of constructing Lyapunov functionals, where H is given in Equation (14), and One can determine the global dynamics of single-group SIR model (21). Applying Theorems 3.1 and 3.2, we obtain the following global dynamics but omit the proof.

Remark 1
The results in Theorem 4.1 improve the corresponding results in [26], which gives part of the proof for this problem when f (S, I) = f (S)I in system (9).

Numerical simulation
In this section, we carry out numerical simulation to illustrate and support our analytical results.
We consider a simpler case in which all groups share the same natural death rate: d i = d for i = 1, 2, . . . , n. Further, we denote d t P i (t − ξ) = −g(t − ξ) and assume that the functions g i (ξ ) are disease-specific only, implying that h i (ξ ) = g(ξ ) for i = 1, 2, . . . , n. We choose the gamma distribution: where b > 0 is a real number and n > 1 is an integer. This is widely used and can approximate several frequently used distributions. For example, when b → 0 + , h n,b (s) will approach the Dirac delta function, and when n = 1, h n,b (s) is an exponentially decaying function. Following the technique and method in [35], defineb which can absorb the exponential term e −δu into the delay kernel. The second equation in Equation (5) can be rewritten as For l = 1, . . . , n, let Thus, for l ∈ {2, . . . , n}, we obtaiṅ For l = 1, we have It follows thatẏ Thus, the integro-differential system (5) is equivalent to the ordinary differential equations , . . .

Summary and discussion
We present a complete mathematical analysis for global asymptotic stability of unique equilibrium P * of system (9); this complements our earlier work [34], where we assumed f ij (S i , I j ) = f i (S i )I j . Theorems 3.1 and 3.2 are also applicable to system (5) with other special nonlinear incidence rates appearing in the literature including f ij (S i , I j ) = S i I j [12,13,19]; f ij (S i , I j ) = ηS i I j /(1 + θS i ) [1,11]; and f ij (S i , I j ) = S q i I j [24,25,37]. The new model includes many existing ones as special cases. For n = 1, system (5) reduces to a single-group SIR model with general relapse distribution and nonlinear incidence rate. Biologically, Theorems 3.1 and 3.2 imply that, if the basic reproduction number R 0 ≤ 1, then the disease  (25) with R 0 = 0.144355 < 1, hence P 0 = (3, 0, 0, 0, 3, 0, 0, 0) is globally stable. Graphs (a) and (b) illustrate that S 1 (t) and S 2 (t) will eventually lead to steady state. Graphs (c) and (d) illustrate that I 1 (t) and I 2 (t) will eventually lead to zero. Initial value is S 1 (0) = 6, S 2 (0) = 4, y 1,1 (0) = 0.1, y 1,2 (0) = 0.1, y 2,1 (0) = 0.1, y 2,2 (0) = 0.1, I 1 (0) = 2, I 2 (0) = 1.  (25) with R 0 = 1.13422 > 1, hence P * = (0.621928, 3.82371, 0.377123, 0.377123, 0.621928, 3.82371, 0.377123, 0.377123) is globally stable. Graphs (a) and (b) illustrate that S 1 (t) and S 2 (t) will eventually lead to steady state. Graphs (c) and (d) illustrate that I 1 (t) and I 2 (t) will eventually lead to steady state. Initial value is S 1 (0) = 6, S 2 (0) = 4, y 1,1 (0) = 0.1, y 1,2 (0) = 0.1, y 2,1 (0) = 0.1, y 2,2 (0) = 0.1, I 1 (0) = 2, I 2 (0) = 1. always dies out from all groups; if R 0 > 1, then the disease always persists in all groups at the unique EE level, irrespective of the initial conditions. On the other hand, Theorems 3.1, 3.2 and 4.1 demonstrate that heterogeneity and nonlinear incidence rate do not alter the dynamical behaviour of the SIR model with general relapse distribution and nonlinear incidence rate. Compared to results in [13,14], the group structure in system (9) greatly increases the complexity exhibited in the derivatives of the Lyapunov functionals. The key to our analysis is a complete description of the patterns exhibited in the derivative of the Lyapunov functionals using graph theory. Our approach may provide a frame work for dynamics of sexually transmitted diseases with relapse distribution and multi-group structure. Heterogeneity in the host population can result from different behaviours (e.g. numbers of sexual partners for some sexually transmitted infections). The global dynamics exclude the existence of Hopf bifurcation leading to sustained oscillatory solutions.
We should point here that this work is motivated by Yuan et al. [35] and Sun et al. [31] where disease with latency spreading in a heterogeneous host population and nonlinear incidence of infection and nonlinear removal functions between compartments was considered. In the proof, we utilize a graph-theoretical approach to the method of global Lyapunov functions that is motivated by the works in [13,14,18,19,22,23,26,31,34].