 Research
 Open Access
 Published:
Modelling infectious diseases with relapse: a case study of HSV2
Theoretical Biology and Medical Modellingvolume 14, Article number: 13 (2017)
Abstract
Background
Herpes Simplex Virus Type 2 (HSV2) is one of the most common sexually transmitted diseases. Although there is still no licensed vaccine for HSV2, a theoretical investigation of the potential effects of a vaccine is considered important and has recently been conducted by several researchers. Although compartmental mathematical models were considered for each special case in the previous studies, as yet there are few global stability results.
Results
In this paper, we formulate a multigroup SVIRI epidemic model for HSV2, which enables us to consider the effects of vaccination, of waning vaccine immunity, and of infection relapse. Since the number of groups is arbitrary, our model can be applied to various structures such as risk, sex, and age group structures. For our model, we define the basic reproduction number ℜ_{0} and prove that if ℜ_{0}≤1, then the diseasefree equilibrium is globally asymptotically stable, whereas if ℜ_{0}>1, then the endemic equilibrium is so. Based on this global stability result, we estimate ℜ_{0} for HSV2 by applying our model to the risk group structure and using US data from 2001 to 2014. Through sensitivity analysis, we find that ℜ_{0} is approximately in the range of 23. Moreover, using the estimated parameters, we discuss the optimal vaccination strategy for the eradication of HSV2.
Conclusions
Through discussion of the optimal vaccination strategy, we come to the following conclusions. (1) Improving vaccine efficacy is more effective than increasing the number of vaccines. (2) Although the transmission risk in female individuals is higher than that in male individuals, distributing the available vaccines almost equally between female and male individuals is more effective than concentrating them within the female population.
Background
Herpes Simplex Virus Type 2 (HSV2) is one of the most common sexually transmitted diseases, and has infected about 417 million people aged 1549 worldwide [1]. Although there is still no licensed vaccine for HSV2, a theoretical investigation of the potential effects of a vaccine is considered important and has recently been conducted by several researchers (see [2–4]). In [2, 3], compartmental epidemic models with vaccination for HSV2 were considered and the effectiveness of the vaccination was discussed in connection with the basic reproduction number ℜ_{0} (see [5]) through numerical simulations. However, there was little discussion about the stability of each equilibrium. As observed in several papers on epidemic models with vaccination (see, for instance, [6–8]), backward bifurcation can occur at ℜ_{0}=1 for some special models and ℜ_{0}<1 does not necessarily imply the global asymptotic stability of the diseasefree equilibrium, that is, the eradication of the disease. In that case, the vaccination effort solely to make ℜ_{0}<1 has less significance. Therefore, a global stability analysis is critical for theoretically justifying the epidemiological discussion.
In [4], Lou et al. considered a compartmental epidemic model for HSV2 with age and risk group structures and discussed the effectiveness of the vaccination together with the global stability analysis of each equilibrium. In their study, the vaccination was limited to female individuals, who are known to be the highrisk group for HSV2, and it was concluded that such a vaccination strategy can reduce the total infections in both females and males. However, to support their conclusion, we need to consider a more general model in which male individuals can also benefit from the vaccination and show that the optimal distribution ratio of the vaccines is 1 to 0 for female and male individuals. In this paper, we consider such a general model and investigate the optimal distribution ratio of the vaccines. As opposed to their conclusion, our result shows that distributing the vaccines almost equally to females and males is more effective for the eradication of HSV2 than concentrating them within the female population.
To consider the effect of vaccination with imperfect immunity, SVIR epidemic models are often formulated, in which the total population is subdivided into the susceptible (S), vaccinated (V), infective (I) and recovered (R) populations (see, for instance, [2, 6–10]). However, to take into account the relapse of HSV2 (see [2, 11]), it is necessary to also consider a direct transition from R to I. Thus, in this paper, we formulate a multigroup SVIRI epidemic model for HSV2, which enables us to consider the effects of vaccination, of waning vaccine immunity, and of infection relapse. Since the number of groups is arbitrary, our model can be applied to various structures such as risk, sex, and age group structures. In the empirical portion of this paper, we apply our model to the risk group structure and estimate the basic reproduction number ℜ_{0} for HSV2 by using data from the US from 2001 to 2014. Since the infective population of HSV2 seems to be in endemic equilibrium, the estimation of ℜ_{0} must be carried out under the global asymptotic stability of the endemic equilibrium. However, in general, the global asymptotic stability of the endemic equilibrium is not trivial.
Recently, multigroup epidemic models have been studied by many authors [10, 12–24]. One of the most effective approaches for global stability analysis of multigroup epidemic models is the graphtheoretic approach developed by Guo et al. [14]. Since our model has a quite complex form with the paths from V to S (the waning of vaccineinduced immunity), R to I (relapse) and distributed time delay, the global asymptotic stability analysis is challenging. In this paper, by applying the graphtheoretic approach as in [14] together with an approach of max function as in [10], we prove that if ℜ_{0}>1, then the endemic equilibrium is globally asymptotically stable, whereas if ℜ_{0}≤1, then the diseasefree equilibrium is so. Based on this theoretical result, we estimate ℜ_{0} for HSV2 by using US data from 2001 to 2014. By using the estimated parameters, we discuss the optimal vaccination strategy for the eradication of HSV2.
Methods
The general multigroup SVIRI epidemic model
Let \(n \in \mathbb {N}\) be the number of groups and let \({\mathcal {N}} := \left \{ 1,2,\cdots, n \right \}\). Let N _{ i }(t) be the sexually active population in group \(i \in {\mathcal {N}}\) at time t. Let us divide N _{ i }(t) into four subpopulations: susceptible S _{ i }(t), vaccinated V _{ i }(t), infective I _{ i }(t), and recovered R _{ i }(t). Thus, N _{ i }(t)=S _{ i }(t)+V _{ i }(t)+I _{ i }(t)+R _{ i }(t) for all \(i \in {\mathcal {N}}\). We make the following assumptions: (A1) The number of individuals becoming sexually active in group \(i \in {\mathcal {N}}\) per unit time is b _{ i }>0. (A2) The per capita rate of removal from the sexual activity in group \(i \in {\mathcal {N}}\) is μ _{ i }>0. (A3) The coefficient for disease transmission from infective individuals in group \(j \in {\mathcal {N}}\) to uninfected (susceptible or vaccinated) individuals in group \(i \in {\mathcal {N}}\) is β _{ ij }≥0. The matrix \((\beta _{ij})_{i,j \in {\mathcal {N}}}\) is irreducible. The vaccine efficacy in group \(i \in {\mathcal {N}}\) is σ _{ i }∈ [ 0,1] and the force of infection to vaccinated individuals in group \(i \in {\mathcal {N}}\) is weakened by multiplying σ _{ i }. That is,
are the forces of infection to susceptible and vaccinated individuals in group \(i \in {\mathcal {N}}\) at time t≥0, respectively. Here we assume standard incidence. (A4) The per capita vaccination rate for susceptible individuals in group \(i\in {\mathcal {N}}\) is v _{ i }>0. The per capita rate for the waning of vaccineinduced immunity for vaccinated individuals in group \(i \in {\mathcal {N}}\) is ω _{ i }≥0. (A5) The per capita recovery rate of infective individuals in group \(i \in {\mathcal {N}}\) is γ _{ i }>0. (A6) The survival probability for recovered individuals in group \(i \in {\mathcal {N}}\), who spent time t in the recovered class, is \(P_{i}(t):=\exp (\int _{0}^{t} \delta _{i}(\eta) \mathrm {d}\eta)\), where δ _{ i }(η) denotes the relapse risk for individuals who spent time η in the recovered class in group i. For each \(i \in {\mathcal {N}}\), \(\delta _{i} \in L_{\text {loc}, +}^{1} (0,+\infty)\) and \(\int _{0}^{+\infty }\delta _{i}(\eta)\mathrm {d}\eta = +\infty \).
Under assumptions (A1)(A2), we see that the time variation of \(N_{i}(t), \ i \in {\mathcal {N}}\) is governed by the following differential equation:
From the variation of constants formula, we easily see that \({\lim }_{t\to +\infty } N_{i}(t) = b_{i}/\mu _{i} =: N_{i}^{*}, \ i \in {\mathcal {N}}\). Hence, without loss of generality, we can assume that \(N_{i}(t) \equiv N_{i}^{*}, \ i \in {\mathcal {N}}\). Then, under assumptions (A1)(A4), we obtain the differential equations for S _{ i }(t) and V _{ i }(t), \(i \in {\mathcal {N}}\) as follows:
Under assumptions (A5)(A6), the recovered population in group \(i \in {\mathcal {N}}\) at time t is given by
By differentiating (4), we obtain the following integrodifferential equation for R _{ i }(t), \(i \in {\mathcal {N}}\).
From (1)(5) we obtain the integrodifferential equation for I _{ i }(t), \(i \in {\mathcal {N}}\) as follows.
Under this setting, we arrive at the following main model in this paper.
Note that the differential equation of \(R_{i}(t), \ i \in {\mathcal {N}}\) can be omitted since it does not appear in the above three equations.
The equilibria of system (6) can be obtained as the solution of the following algebraic equations.
where
Note that
Hence, we have μ _{ i }+γ _{ i }−Q _{ i }>0 for all \(i \in {\mathcal {N}}\).
It is easy to see that the trivial solution of (7) such that I _{ i }=0 for all \(i \in {\mathcal {N}}\) always exists. It is called the diseasefree equilibrium of system (6) and we write it as \(E^{0} := \left (S_{1}^{0}, V_{1}^{0}, 0, S_{2}^{0}, V_{2}^{0}, 0, \cdots, S_{n}^{0}, V_{n}^{0}, 0 \right) \in \mathbb {R}_{+}^{3n}\), where
Existence of the endemic equilibrium \(E^{*} : = \left (S_{1}^{*}, V_{1}^{*}, I_{1}^{*}, \cdots, S_{n}^{*}, V_{n}^{*}, I_{n}^{*} \right) \in \mathbb {R}_{+}^{3n}\) such that \(I_{i}^{*} > 0\) for all \(i \in {\mathcal {N}}\) will be discussed in connection with the basic reproduction number ℜ_{0}, which is defined as the expected number of secondary cases produced by a typical infected individual during its entire period of infectiousness at the initial invasion phase into a fully susceptible population, and given by the spectral radius of the next generation matrix (see [25]). Let
Then, according to [25], the next generation matrix is given by
Hence, the basic reproduction number ℜ_{0} is defined by
where ρ(·) denotes the spectral radius of a matrix. We will obtain the global stability results for (6) in connection with ℜ_{0} (see the “Results” section).
The special multigroup SVIRI epidemic model for HSV2
The general model (6) can be applied to analyze the field data of HSV2 epidemics. Similar to other sexually transmitted infections, the risk factor for HSV2 infection is sexual behavior. To describe the heterogeneity of HSV2 infection risk between host individuals, we characterize the group as the combination of sex and their sexual behavior. We consider the following levels of sexual activity: x=0,1,2,⋯,5 meaning the number of sexual partners within the last 12 months, where x=5 implies the number of sexual partners is 5 or more. Let y∈{1,2} denote the sex, 1 denotes male and 2 denotes female. Then, the risk group is characterized by i∈{1,2,⋯,12} in the following way: i=2x _{ i }+y _{ i }, where
For example, i=2 corresponds to (x _{ i },y _{ i })=(0,2) and implies the group of female individuals with no sexual partners and i=11 corresponds to (x _{ i },y _{ i })=(5,1) and implies the group of male individuals with 5 or more sexual partners. Then, (6) can be written as follows:
Note that (9) is a special case of (6). In this section, we assume that δ _{ i }(ξ)≡δ _{ i }>0 for all i∈{1,2,⋯,12}. Note that the assumption (A6) is satisfied. In this case, we have:
Hence, together with the Eq. 5 of R _{ i }(t), (9) can be simplified to the following multigroup SVIRI epidemic model.
No vaccine against HSV2 infection is currently available, so we ignore the vaccinated class V _{ i }, i∈{1,2,⋯,12} in the estimation of ℜ_{0}. Then, (10) can be rewritten as follows.
The basic reproduction number ℜ_{0} for (11) is obtained as the spectral radius of the following matrix.
where Q _{ i }=δ _{ i } γ _{ i }/(μ _{ i }+δ _{ i }) and we write β _{ ij } as β _{ i,j } for improved readability.
Transmission rates between the risk groups i and j, β _{ ij }, depend on sexual behavior and sex. We modeled β _{ ij } as follows;
The meaning of each symbol for β _{ ij } is as follows.

\(\rho _{x_{i} y_{i}}\) denotes the HSV2 infection risk for the risk group i. The risk group is stratified by sex and the number of partners within the last 12 months, the risk group i denotes the individuals whose number of partners within the last 12 months is x _{ i } and sex is y _{ i }. \(\rho _{x_{i} y_{i}}\) is given by;
$$\rho_{x_{i} y_{i}} = c_{y_{i}} (x_{i} + 1)^{\phi}. $$Here, similar to previous modelling studies of sexually transmitted infections, we modeled the relationship between infection risk and sexual behavior by a power law function [26].

c denotes the sex specific HSV2 transmission coefficient.

ϕ denotes the exponent parameter describing the heterogeneity of the infection risk between different sexual behaviors.

R denotes the mixing matrix between the risk groups defined by sexual behavior, x;
$$\mathbf{R}_{x_{i} x_{j}} = q \delta_{x_{i} x_{j}} + (1q) \frac{\sum_{y} \rho_{xy} N_{xy}}{\sum_{x} \sum_{y} \rho_{xy} N_{xy}}. $$This is the classical oneparameter ‘preferred mixing’ formulation, proposed by [27].

δ denotes Kronecker’s delta.

q denotes assortative coefficient. When q=0, the mixing between risk groups defined by sexual behavior is ‘proportionately mixing’, and the mixing is ‘fully assortative mixing’ when q=1.

S denotes the mixing matrix between sexes;
$$\mathbf{S} = \left(\begin{array}{cc} a & 1 a \\ 1a & a \end{array} \right). $$ 
a denotes the proportion of homosexual behavior.
We will use the special model (11) with transmission rate (12) to estimate the basic reproduction number ℜ_{0} for HSV2 (see the “Results” section), and (10) with (12) to discuss the effectiveness of vaccination strategy (see the “Discussion” section).
Results
The main theorem
The main theorem of this paper is obtained for the general multigroup SVIRI epidemic model (6). Since (6) has an infinite time delay, we define the fading memory space (see, for instance, [28, 29]) as follows:
where Δ is a positive constant such that \(0 < \Delta < \min _{i \in {\mathcal {N}}} \{ \mu _{i} \}\). Let us define the following state space for system (6):
The following proposition is proved:
Proposition 1
Ω is positively invariant for system (6).
The main theorem of this paper is as follows.
Theorem 1
Let ℜ_{0} and Ω be defined by (8) and (14), respectively. Let \(\bar {\Omega }\) denote the closure of Ω. (i) If ℜ_{0}≤1, then the diseasefree equilibrium \(E^{0} \in \bar {\Omega }\) of system (6) is globally asymptotically stable in Ω and there exists no endemic equilibrium E ^{∗} in \(\bar {\Omega }\). (ii) If ℜ_{0}>1, then the system (6) has the unique endemic equilibrium E ^{∗} in Ω and it is globally asymptotically stable in Ω.
For the proofs of Proposition 1 and Theorem 1, see the Appendix.
Theorem 1 still works for (10) since it is a special case of (6). In particular, although (10) does not include the integrated time delay, to our knowledge, there is no previous study on the global asymptotic stability of the endemic equilibrium of model (10). From this viewpoint, our main theorem can be regarded as valuable for the empirical study in the subsequent sections.
Estimation of ℜ_{0} for HSV2
Based on Theorem 1, we estimate the basic reproduction number ℜ_{0} for HSV2 in the US from 2001 to 2014. For the estimation of ℜ_{0}, we use the special model (11) with transmission rate (12). Note that (11) corresponds to the case where v _{ i }=σ _{ i }=ω _{ i }=0 for all i∈{1,2,⋯,12}. Although the case where v _{ i }=0 for all i∈{1,2,⋯,12} is excluded under assumption (A4), it is easy to check in a completely similar way as in the Appendix that the global stability result similar to Theorem 1 holds.
Previous study derived the value of δ _{ i } and γ _{ i } from empirical data, δ _{ i } and γ _{ i } are parameterized based on [30], 1/δ _{ i }=78.5 days and 1/γ _{ i }=13 days for all i∈{1,2,⋯,12}. Here note that we can regard μ _{ i } as the removal rate from our system, which is given by the sum of the sexualinactivation rate and the mortality rate among those who are sexually active. We assume that the sexual life span is 50 years (1565 years old) and parameterize the mortality rate by the national representative census data in the US [31], μ _{ i }=0.0231 per year for all i∈{1,2,⋯,12}. Based on the previous studies [32] and [33], we obtain the estimations q=0.3 and a=0.02, respectively (see Table 1).
Using the observed data of seroprevalence of HSV2 in the US from 2001 to 2014 reported by [34], sex specific transmission coefficient c and the exponent parameter ϕ were estimated by maximum likelihood estimation. Since the antibody against HSV2 infection (IgG) provides lifelong immunity [35], we fitted I+R to the observed data of the number of seropositive cases for the estimation of c and ϕ. To estimate c and ϕ, endemic equilibria of I _{ i } and R _{ i } were solved numerically with varied c _{1} and c _{2} and ϕ, and the set of c _{1} and c _{2} maximizing the likelihood function was explored. The likelihood function for c _{1} and c _{2} is given by
Here pmf denotes the probability mass function, bin denotes a binomial distribution, \(N^{\text {data}}_{i,T}\) denotes the observed data of the size of the risk group i in sampling year T, and \(P^{\text {data}}_{i,T}\) denotes the observed data of the number of HSV2seropositive cases in the risk group i in sampling year T, respectively. For the confidence interval (CI) of the estimated parameter, profile likelihoodbased confidence intervals were calculated. Using estimated c and ϕ the basic reproduction number ℜ_{0} for HSV2 in the US was calculated. Figure 1 shows the comparison between the observed data of seroprevalence of HSV2 and the model estimates. The estimated c are, transmission coefficient for male, c _{1}=0.228 (95% CI 0.225 to 0.231), transmission coefficient for female, c _{2}=1.78 (95% CI 1.75 to 1.81), exponent parameter ϕ=0.700 (95% CI 0.693 to 0.707) and estimated ℜ_{0}=2.07 (95% CI 2.03 to 2.11).
Sexual behavior shows wide variation between host individuals. To assess the sensitivity of sexual behavior to ℜ_{0} of HSV2, we conducted a sensitivity analysis of the parameters describing sexual behavior, i.e., the proportion of homosexual partnership a and assortativity coefficient for the mixing between risk groups q. Fig. 2 shows the relation of a and q to estimated ℜ_{0}, ℜ_{0} increase if i) a increases, and ii) q decreases. The realistic variations of a and q [36, 37] can induce the variation of ℜ_{0}, which is approximately demonstrated in the range of 23.
Discussion
Using the demographic and epidemiological parameters obtained above, we discuss the effectiveness of each vaccination strategy. We investigate the sensitivity of the basic reproduction number ℜ_{0} to the vaccination parameters, that is, the vaccination rate among susceptible population v and the vaccination efficacy σ. Here we have assumed that vaccination is conducted with the same rate v for the susceptible population over time. For simplicity, we assume that the efficacy of vaccine σ is the same for all risk groups.
We first consider the case that vaccination rate v is the same between males and females. In this case, the basic reproduction number ℜ_{0} with different σ when v varies over (0,1) is shown in Fig. 3. We see from Fig. 3 that, if σ is 0.3 or smaller, ℜ_{0} can be less than 1. On the other hand, if σ is 0.4 or larger, ℜ_{0} cannot be less than 1 for any v∈(0,1). This implies that decreasing σ is more important than increasing v to reduce the basic reproduction number ℜ_{0}. That is, improving the vaccine efficacy is more important for the eradication of HSV2 than increasing the number of vaccines.
We next discuss the optimal sex ratio of the vaccinated population to control HSV2. HSV2 infection is observed among females more frequently than males, “opportunistic” vaccination can induce higher vaccination coverage among females than males. To assess the optimal sex ratio of the vaccination rate, we expand the vaccination rate v as follows;
Here p denotes the sex ratio of vaccination. Figure 4 shows the relationship between p, σ and ℜ_{0}, we assumed v=0.9 as the representative value. Interestingly, small or large p increases ℜ_{0}. This implies that vaccination biased to females (small p) or males (large p) can result in persistence of the disease. In particular, it is noteworthy that the curves in Fig. 4 are almost symmetric with respect to p and the minimum is attained near the center p=0.5. This implies that vaccination distributed equally to females and males is optimal for the eradication of the disease even though the transmission coefficient for males is lower than that for females.
Conclusion
In this paper, we have formulated the multigroup SVIRI epidemic model (6), which enables us to consider the effects of vaccination, the waning of vaccineinduced immunity, and relapse. We have defined the basic reproduction number ℜ_{0} and proved Theorem 1, which states that if ℜ_{0}≤1, then the diseasefree equilibrium E ^{0} is globally asymptotically stable, whereas if ℜ_{0}>1, then the endemic equilibrium E ^{∗} is so. Based on Theorem 1, we have estimated the basic reproduction number ℜ_{0} for HSV2 as 2.07 (95% CI 2.03 to 2.11) by using US HSV2 data from 2001 to 2014. Through the sensitivity analysis for uncertain parameters on sexual behavior, we have found that ℜ_{0} is approximately in the range of 23. Furthermore, using sensitivity analysis for vaccination parameters, we have discussed the effectiveness of the vaccination. As a result, we have come to the following conclusions. (1) Improving vaccine efficacy is more effective than increasing the number of vaccines. (2) Although the transmission risk in female individuals is higher than that in male individuals, distributing vaccines almost equally to females and males is more effective than concentrating them within the female population.
Appendix
Proof of Proposition 1
We first show the positivity of the solution of system (6). Suppose that there exist t _{1}>0 and \(i^{*} \in {\mathcal {N}}\) such that S _{ i }(t)>0 and V _{ i }(t)>0 for all t∈[0,t _{1}) and \(i \in {\mathcal {N}}\) and \(\min \left (S_{i^{*}}(t_{1}), V_{i^{*}}(t_{1}) \right) =0\). By the variation of constants formula, we have from the first equation in the system (6) that
Hence, \(V_{i^{*}}(t_{1}) = 0\phantom {\dot {i}\!}\). However, by the variation of constants formula, we have from the second equation in the system (6) that
which is a contradiction. Hence, we see that S _{ i }(t)>0 and V _{ i }(t)>0 for all t>0 and \(i \in {\mathcal {N}}\).
Suppose that there exist t _{2}>0 and \(\tilde {i} \in {\mathcal {N}}\) such that I _{ i }(t)>0 for all t∈[0,t _{2}) and \(i \in {\mathcal {N}}\) and \(I_{\tilde {i}} (t_{2}) = 0\). By the variation of constants formula, we have from the third equation in the system (6) that
where \(h_{i} (t) := \int _{0}^{+\infty } \delta _{i}(\xi) \gamma _{i} I_{i}(t\xi)e^{\mu _{i} \xi } e^{\int _{0}^{\xi } \delta _{i}(\eta) \mathrm {d}\eta } \mathrm {d}\xi \). We see that h _{ i }(t)≥0 for all \(i \in {\mathcal {N}}\) and t∈ [ 0,t _{1}). Hence, from (15), we have \(I_{\tilde {i}}(t_{2}) > 0\), which is a contradiction. Hence, we see that I _{ i }(t)>0 for all t>0 and \(i \in {\mathcal {N}}\).
The boundedness of the solution of system (6) immediately follows from the fact that N i′(t)=b _{ i }−μ _{ i } N _{ i }(t), S i′(t)≤b _{ i }−(μ _{ i }+v _{ i })S _{ i }(t)+ω _{ i } V _{ i }(t) and V i′(t)≤v _{ i } S _{ i }(t)−(μ _{ i }+ω _{ i })V _{ i }(t) for all t>0 and \(i \in {\mathcal {N}}\). This completes the proof.
Proof of (i) of Theorem 1
We define the following matrix, which corresponds to the next generation matrix:
In fact, it is easy to see that \(\rho (M^{0}) = \rho (\mathcal {K}) =\Re _{0}\).
First we show that system (6) has no endemic equilibrium \(E^{*} \in \bar {\Omega }\). Let us define the following matrixvalued function on \(\mathbb {R}^{2n}\), which is equal to M ^{0} if \((S_{1},V_{1},\cdots,S_{n},V_{n}) = (S^{0}_{1},V^{0}_{1}, \cdots S^{0}_{n}, V^{0}_{n})\):
Suppose that \((S_{1},\cdots,S_{n}) \neq (S_{1}^{0},\cdots, S_{n}^{0})\). Then, from assumptions (A1)(A6), we see that 0<M(S _{1},V _{1},⋯,S _{ n },V _{ n })<M ^{0}, where 0 denotes the zero matrix and the inequality implies that it holds for each element and each of the two matrices are not equal. Then, since it follows from assumptions (A1)(A6) that matrices M ^{0} and M ^{0}+M(S _{1},V _{1},⋯,S _{ n },V _{ n }) are nonnegative and irreducible, we can apply the PerronFrobenius theorem (see [38, Corollary 2.1.5]) to obtain that ρ(M(S _{1},V _{1},⋯,S _{ n },V _{ n }))<ρ(M ^{0})≤1. This implies that the equation M(S _{1},V _{1},⋯,S _{ n },V _{ n })(I _{1},⋯I _{ n })^{T}=(I _{1},⋯I _{ n })^{T} has only the trivial solution (I _{1},⋯,I _{ n })^{T}=0, where T denotes the transpose of a vector. This implies that E ^{∗} does not exist in \(\bar {\Omega }\).
Next we show the global asymptotic stability of E ^{0}. It follows from the PerronFrobenius theorem (see [38, Theorem 2.1.4]) that M ^{0} has a strictly positive left eigenvector (ℓ _{1},⋯,ℓ _{ n }), ℓ _{ i }>0, \(i \in {\mathcal {N}}\) corresponding to the eigenvalue ρ(M ^{0}): (ℓ _{1},⋯,ℓ _{ n })ρ(M ^{0})=(ℓ _{1},⋯,ℓ _{ n })M ^{0}. Let \(c_{i} : = \ell _{i} /\left (\mu _{i} + \gamma _{i}  Q_{i} \right), \ i \in {\mathcal {N}}\) and \(J_{i} (t) : = \int _{t}^{+\infty } \delta _{i}(\xi) \gamma _{i} e^{\mu _{i} \xi } e^{\int _{0}^{\xi } \delta _{i}(\eta) \mathrm {d}\eta } \mathrm {d}\xi, \ i \in {\mathcal {N}}\) and consider the following Lyapunov function.
From assumption (A6), J _{ i }(t)≥0, \(i \in {\mathcal {N}}\) for all t≥0 and hence, \(\mathcal {L}_{DFE} \geq 0\) and the equality holds if and only if (I _{1},⋯,I _{ n })≡0. Note that
Hence, the derivative of \(\mathcal {L}_{DFE}\) gives
where E _{ n } denotes the ndimensional unit matrix and · denotes the product of vectors. It is easy to see that when ℜ_{0}<1, \(\mathcal {L}_{DFE}' = 0\) holds if and only if I _{ i }=0 for all \(i \in {\mathcal {N}}\), that is, the solution is in the diseasefree equilibrium E ^{0}. When ℜ_{0}=1, from the third equality in (17), we see that \(\mathcal {L}_{DFE}' = 0\) implies
Suppose that \((S_{1},V_{1},\cdots,S_{n},V_{n}) \neq \left (S_{1}^{0},V_{1}^{0},\cdots,S_{n}^{0},V_{n}^{0}\right)\). Then (ℓ _{1},⋯,ℓ _{ n })·M(S _{1},V _{1},⋯,S _{ n },V _{ n })<(ℓ _{1},⋯,ℓ _{ n })·M ^{0}=ρ(M ^{0})(ℓ _{1},⋯,ℓ _{ n })=(ℓ _{1},⋯,ℓ _{ n }). Hence, (18) has only the trivial solution such that I _{ i }=0 for all \(i \in {\mathcal {N}}\). This implies that \(\mathcal {L}_{DFE}' = 0\) holds only in the diseasefree equilibrium \(E^{0} \in \bar {\Omega }\). Consequently, from the LaSalle’s invariance principle (see [39]), we can conclude that the diseasefree equilibrium E ^{0} is globally asymptotically stable.
Proof of (ii) of Theorem 1
If ℜ_{0}>1, then \(\left (\ell _{1}, \cdots, \ell _{n} \right) \cdot \left (M \left (S_{1}^{0},V_{1}^{0}, \cdots, S_{n}^{0}, V_{n}^{0} \right)  E_{n} \right) \cdot (I_{1},\cdots, I_{n})^{T} = \left (\rho (M^{0})  1 \right) \left (\ell _{1}, \cdots, \ell _{n} \right) \cdot (I_{1},\cdots, I_{n})^{T} \ > \ 0\). Hence, we see from the third equality in (17) that in a neighborhood of \(\left (S_{1}^{0},V_{1}^{0},\cdots, S_{n}^{0}, V_{n}^{0}\right)\), \(\mathcal {L}_{DFE}' > 0\). This implies the instability of the diseasefree equilibrium E ^{0}.
Since the diseasefree equilibrium of E ^{0} of system (6) is unstable if ℜ_{0}>1, we see from the uniform persistence result of [40] and an argument as in the proof of Proposition 3.3 of [41] that system (6) is uniformly persistent. That is, there exists a positive constant c>0 such that for any initial value, it holds that \(\liminf _{t\to +\infty } S_{i}(t) \geq c, \ \liminf _{t\to +\infty } V_{i}(t) \geq c\) and \(\liminf _{t\to +\infty } I_{i}(t) \geq c\) for all \(i \in {\mathcal {N}}\). Since the uniform persistence together with the uniform boundedness implies the existence of an interior equilibrium (see [42, 43]), we see that system (6) has an endemic equilibrium E ^{∗}∈Ω. From (7), we see that the components \(\left (S_{1}^{*},V_{1}^{*},I_{1}^{*},\cdots, S_{n}^{*}, V_{n}^{*}, I_{n}^{*}\right)\) of E ^{∗} satisfy the following equations.
As in [14], we define \(\theta _{ij} := \left (S_{i}^{*} + \sigma _{i} V_{i}^{*} \right) \beta _{ij} I_{j}^{*}/N_{j}^{*}\), \(i, j \in {\mathcal {N}}\) and
Let φ:=(φ _{1},⋯,φ _{ n })^{T} be a basis of the solution space of linear system Θ φ=0. Then, from [14, Lemma 2.1], we see that the dimension of the solution space is 1 and φ _{ i }>0, \(i \in {\mathcal {N}}\). In particular, from the form of matrix Θ, the following equality holds.
Using this φ and H(x):=x−1− lnx≥H(1)=0, we consider the following Lyapunov functional to prove the global asymptotic stability of E ^{∗}.
In order to make this function welldefined, without loss of generality, we can restrict our attention to the solution such that \(I_{i}(s)=\varphi _{i} (s), \ i \in \mathcal {N}\) on (−∞,0], where φ _{ i }(0)=I _{ i }(0) and \(0 < m_{i} < \varphi _{i}(s) < M_{i} < +\infty, \ s \in (\infty, 0], \ i \in \mathcal {N}\) for positive constants m _{ i } and M _{ i }, \(i \in \mathcal {N}\). Then, from the positive invariance of set Ω and the uniform persistence of system (6), we see that the Lyapunov functional \(\mathcal {L}_{EE}\) is welldefined.
Using (19), we can calculate the derivative of \(\mathcal {L}_{EE}\) as follows.
Now it follows from integration by parts that
Hence, (23) can be calculated as follows:
Hence, using (20), (21) and (25), we can calculate (24) as follows:
Using the inequality of arithmetic and geometric means, we see that the first three terms in the righthand side of (26) are nonpositive and equal to zero if and only if \(\left (S_{i}, V_{i} \right) = \left (S_{i}^{*}, V_{i}^{*} \right)\), \(i \in {\mathcal {N}}\). From the positivity of the function H(x), we see that the last term in the righthand side of (26) is nonpositive. Hence, taking the maximum as in [10] and using the graphtheoretic approach as in [14], we can evaluate (26) as follows:
where Γ denotes the set of all unicyclic graphs included in directed graphs with vertices {1,2,⋯,n}, G denotes the unicyclic graph included in Γ, w(G) denotes the weight of graph G, CG denotes the unicycle included in G and A(C G) denotes the set of all arcs included in CG. For instance, for a unicycle C G:1→2→1, we have A(C G)={(1,2),(2,1)} and thus,
We see that all elements in the max in the last expression of the above formula are non positive because of the inequality of arithmetic and geometric means. Similarly, we can easily check that for all unicycles CG with at most n vertices, the second sum in the last expression of (27) are nonpositive (see [10, Proof of Theorem 4.1]). Hence, \(\mathcal {L}_{EE}'\) is nonpositive and it is easy to check that the equality \(\mathcal {L}_{EE}' = 0\) holds if and only if \(\left (S_{1},V_{1},I_{1},\cdots, S_{n},V_{n},I_{n} \right) = \left (S_{1}^{*},V_{1}^{*},I_{1}^{*},\cdots, S_{n}^{*},V_{n}^{*},I_{n}^{*} \right) \). This implies, from the LaSalle’s invariance principle, that the endemic equilibrium E ^{∗} is globally asymptotically stable.
References
 1
World Health Organization. Media centre: Herpes simlex virus. 2017. Available from: http://www.who.int/. Accessed 25 Jan 2017.
 2
Alsallaq RA, Schiffer JT, Longini IM, Wald A, Corey L, AbuRaddad LJ. Population level impact of an imperfect prophylactic HSV2 vaccine. Sex Transm Dis. 2010; 37:290–7.
 3
Blower S. Modelling the genital herpes epidemic. Herpes. 2004; 3:138A–146A.
 4
Lou Y, Qesmi R, Wang Q, Steben M, Wu J, Hefferman JM. Epidemiological impact of a genital herpes type 2 vaccine for young females. PLoS ONE. 2012; 7:e46027.
 5
Diekmann O, Heesterbeek JAP, Metz JAJ. On the definition and the computation of the basic reproduction ratio Ro in models for infectious diseases in heterogeneous populations. J Math Biol. 1990; 28:365–82.
 6
KribsZaleta CM, VelascoHernández JX. A simple vaccination model with multiple endemic states. Math Biosci. 2000; 164:183–201.
 7
Arino J, McCluskey CC, van den Driessche P. Global results for an epidemic model with vaccination that exhbits backward bifurcation. SIAM J Appl Math. 2003; 64:260–76.
 8
Alexander ME, Bowman C, Moghadas SM, Summers R, Gumel AB, Sahai BM. A vaccination model for transmission dynamics of influenza. SIAM J Appl Dynam Syst. 2004; 4:503–24.
 9
Liu X, Takeuchi Y, Iwami S. SVIR epidemic model with vaccination strategies. J Theoret Biol. 2008; 253:1–11.
 10
Kuniya T. Global stability of a multigroup SVIR epidemic model. Nonlinear Anal RWA. 2013; 14:1135–43.
 11
van den Driessche P, Zou X. Modeling relapse in infectious diseases. Math Biosci. 2007; 207:89–103.
 12
Lajmanovich A, Yorke JA. A deterministic model for gonorrhea in a nonhomogeneous population. Math Biosci. 1976; 28:221–36.
 13
Feng Z, Huang W, CastilloChavez C. Global behavior of a multigroup SIS epidemic model with age structure. J Diff Equat. 2005; 218:292–324.
 14
Guo H, Li MY, Shuai Z. Global stability of the endemic equilibrium of multigroup SIR epidemic models. Canada Appl Math Quart. 2006; 14:259–84.
 15
Li MY, Shuai Z, Wang C. Global stability of multigroup epidemic models with distributed delays. J Math Anal Appl. 2010; 361:38–47.
 16
Ding D, Ding X. Global stability of multigroup vaccination epidemic models with delays. Nonlinear Anal RWA. 2011; 12:1991–7.
 17
Kuniya T. Global stability analysis with a discretization appraoch for an agestructured multigroup SIR epidemic model. Nonlinear Anal RWA. 2011; 12:2640–55.
 18
Shu H, Fan D, Wei J. Global stability of multigroup SEIR epidemic models with distributed delays and nonlinear transmission. Nonlinear Anal RWA. 2012; 13:1581–92.
 19
Yuan C, Jiang D, O’Regan D, Agarwal RP. Stochastically asymptotically stability of the multigroup SEIR and SIR models with random perturbation. Commun Nonlinear Sci Numer Simulat. 2012; 17:2501–16.
 20
Muroya Y, Enatsu Y, Kuniya T. Global stability for a multigroup SIRS epidemic model with varying population sizes. Nonlinear Anal RWA. 2013; 14:1693–704.
 21
Wang Z, Fan X, Han Q. Global stability of deterministic and stochastic multigroup SEIQR models in computer network. Appl Math Model. 2013; 37:8673–86.
 22
Zhang L, Pang J, Wang J. Stability analysis of a multigroup epidemic model with general exposed distribution and nonlinear incidence rates. Abstr Appl Math. 2013. 2013(Article ID.354287).
 23
Wang J, Liu X, Pang J, Hou D. Global dynamics of a multigroup epidemic model with general exposed distribution and relapse. Osaka J Math. 2015; 52:117–38.
 24
Wang J, Pang J, Liu X. Modelling diseases with relapse and nonlinear incidence of infection: a multigroup epidemic model. J Biol Dynam. 2015; 8:99–116.
 25
van den Driessche P, Watmough J. Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Math Biosci. 2002; 180:29–48.
 26
Liljeros F, Edling CR, Amaral LAN, Stanly HE, Aberg Y. The web of human sexual contacts. Nature. 2001; 411:907–8.
 27
Jacquez JA, Simon CP, Koopman J, Sattenspiel L, Perry T. Modelling and analyzing HIV transmission: the effect of contact patterns. Math Biosci. 1988; 92:119–99.
 28
Hino Y, Murakami S, Naito T. Functional differential equations with infinite delay. vol. 1473 of Lecture Notes in Mathematics. Berlin: SpringerVerlag Berlin Heidelberg; 1991.
 29
Röst G, Wu J. SEIR epidemiological model with varying infectivity and infinite delay. Math Biosci Eng. 2008; 5:389–402.
 30
AbuRaddad LJ, Schiffer JT, Ashley R, Mumtaz G, Alsallaw RA, Akala FA, et al. HSV2 serology can be predictive of HIV epidemic potential and hidden sexual risk behavior in the Middle East and North Africa. Epidemics. 2010; 2:173–82.
 31
Kochanek KD, Murphy SL, Xu J, TejadaVera B. Deaths: Final data for 2014. Natl Vital Stat Rep. 2016; 65:1–122.
 32
AbuRaddad LJ, Jr IML. No HIV stage is dominant in driving the HIV epidemic in subSaharan Africa. AIDS. 2008; 22:1055–61.
 33
Mumtaz G, Hilmi N, McFarland W, Kaplan RL, Akala FA, Semini I, et al. Are HIV epidemics among men who have sex with men emerging in the Middle East and North Africa?: a systematic review and data synthesis. PLoS Med. 2011; 8:e1000444.
 34
The National Health and Nutrition Examination Survey (NHANES). 2016. Available from:http://www.cdc.gov/nchs/nhanes/. Accessed 25 July 2016.
 35
van Wagoner NJ, III EWH. Herpes diagnostic tests and their use. Curr Infect Dis Rep. 2012; 14:175–84.
 36
Garnett GP, Anderson RM. Contact tracing and the estimation of sexual mixing patterns: the epidemiology of gonococcal infections. Sex Trans Dis. 1993; 20:181–91.
 37
Sell RL, Wells JA, Wypij D. The prevalence of homosexual behavior and attraction in the United States, the United Kingdom and France: Results of national populationbased samples. Arch Sex Behav. 1995; 24:235–48.
 38
Berman A, Plemmons RJ. Nonnegative matrices in the mathematical sciences. New York: Academic Press; 1979.
 39
LaSalle JP. The Stability of Dynamical Systems. Philadelphia: SIAM; 1976.
 40
Freedman HI, Ruan S, Tang M. Uniform persistence and flows near a closed positively invariant set. J Dynam Diff Equat. 1994; 6:583–600.
 41
Li MY, Graef JR, Wang L, Karsai J. Global dynamics of a SEIR model with varying total population size. Math Biosci. 1999; 160:191–213.
 42
Bhatia NP, Szego GP. Dynamical systems: stability theory and applications. SpringerVerlag; 1967.
 43
Smith HL, Waltman P. The Theory of the Chemostat: Dynamics of Microbial Competition. Cambridge: Cambridge University Press; 1995.
Acknowledgements
We would like to thank the editor and anonymous reviewers for their helpful comments to the earlier version of this paper. We would like to thank Dr. Akihiro Ishii for helpful discussions regarding the serological test for HSV2.
Funding
JW was supported by National Natural Science Foundation of China (nos. 11401182, 11471089), Science and Technology Innovation Team in Higher Education Institutions of Heilongjiang Province (No. 2014TD005). HT was supported by GrantinAid for JSPS Fellows from the Ministry of Education, Culture, Sports, Science, and Technology in Japan. TK was supported by GrantinAid for Young Scientists (B) of Japan Society for the Promotion of Science (No. 15K17585) and the program of the Japan Initiative for Global Research Network on Infectious Diseases (JGRID); from Japan Agency for Medical Research and Development, AMED. RO was supported by GrantinAid for Young Scientists (B) of Japan Society for the Promotion of Science (No. 15K19217) and Precursory Research for Embryonic Science and Technology (PRESTO) grant number JPMJPR15E1 from Japan Science and Technology Agency (JST). The authors were supported by JSPS Bilateral Joint Research Project (Open Partnership).
Availability of data and materials
The data that support the findings of this study are available in the National Health and Nutrition Examination Survey, “http://www.cdc.gov/nchs/nhanes/”.
Author information
Affiliations
Contributions
JW and XY formulated the model. HT improved the whole manuscript. TK carried out the theoretical analysis of the model. RO carried out the epidemiological study including the estimation of ℜ_{0} for HSV2. All authors read and approved the final manuscript.
Corresponding author
Correspondence to Toshikazu Kuniya.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Multigroup SVIRI epidemic model
 Relapse
 Basic reproduction number
 Global asymptotic stability
 Herpes Simplex Virus Type 2
 Vaccination