Mathematical analysis of an age-structured HSV-2 model

Abstract

The stability analysis of an age-structured Herpes Simplex Virus type two is investigated. The classes of latent E are considered in the mathematical model. The threshold conditions for stability have been obtained. The model has both locally and globally asymptotically stable disease-free equilibrium point. Moreover, it is shown that the transformed model has a unique endemic equilibrium point whenever the reproduction number $\mathcal{R}_{0}>1$ . The endemic equilibrium of this model is shown to be locally-asymptotically stable. For the case where the effective contact rate and natural death rate are constants, it is shown that the model has a unique endemic equilibrium which is globally-asymptotically whenever the associated reproduction number exceeds unity. Furthermore, it is shown that adding the age-factor to the corresponding autonomous HSV2 model (considered in Podder and Gumel [23] for the case when quiescent infectious individuals can transmit infection at the same rate as the non-quiescent infectious individuals) does not alter the qualitative dynamics of the autonomous system (with respect to the elimination or persistence of the disease).

Keywords

Age-structure HSV-2 equilibria stability Mathematics Subject Classification (2010): Primary: 92B05 92C60;Secondary: 92D30.

1. Introduction

Human life is affected by numerous types of diseases. Some of them display minor symptoms whereas some are horrible. Since the middle of the twentieth century, vaccines and antibiotics have improved the treatment process of the infectious diseases. But it has been discovered that infectious diseases have not been terminated and they become one of the most pertinent reasons for human suffering and mortality [19]. In the past few decades, the understanding of disease spread has increased. Several mathematical models have been studied to quantify and analyze the disease spread and to control their infections [16, 17, 26, 27].

Herpes Simplex is one of the diseases that is caused by viruses (HSV-1 and HSV-2). There are two types of Herpes Simplex Virus (HSV) depending on which part of the body is infected. HSV type one (HSV-1) causes oral infections, which include infection in the face and mouth. The HSV type two (HSV-2) causes genital infections [1]. HSV-2 is classified as a sexually transmitted disease (STD). It spreads worldwide among human populations and becomes a severe public health hazard [38]. The spread of HSV-2 varies in different parts of the world. Its spread is more in the developing countries, such as some African and Asian countries, in comparison to the countries in America and Europe. Also the HSV-2 prevalence is higher in adults, which is above 50% in some countries. Young women are more susceptible to this infection than men. One of the studies in Kisumu, Kenya reported that about 40% of young women (aged 15–19) are infected with this virus [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37].

The importance of considering the age factor in the mathematical model comes from several studies which have been made in the past few years. These studies indicate that age is a major factor in HSV-2 prevalence. Chayavichitsilp et. al. noted that more than 60% of adults are infected by either HVS-1 or HSV-2 [3]. The highest percentage of HSV-2 infections worldwide is to be found among young ladies.

Numerous studies considered age-structured models [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]. Li et al. in 2001 [20] studied the threshold conditions and stability of the age-structured SEIR epidemic model in 2008. Li and Liu [21] investigated the stability of the age-structured epidemiological model for hepatitis C. Safi et al. [25] studied the age-structured SEIR epidemic model along with its prescribed treatment. They extended the work done by Elvebac et al. [4] by including the classes for latent E into the study. This paper includes the age factor in to the mathematical model of HSV-2. This inclusion extends the work done by Podder and Gumel in 2010 [23]. Furthermore, this model extends the HSV-2 models studied by Blower et al. [2] and Gershengorn and Blower [6] by considering the ones who have been exposed to HSV-2 but show no clinical symptoms of the disease classes (E).

This paper is organized as follows: in Section 2 the mathematical mode with some basic properties are introduced, the local and global stabilities of the disease-free equilibrium are investigated in Section 3, the existence and uniqueness of the endemic equilibria of the initial value problem is studied and the local stability of the endemic equilibrium point is explored in Section 4.

2. Model formulation

The total population, at time $t$ and age $a$ (where 0 $<a\leqslant A$ , with $A$ being the maximum possible age attained by individuals in the population) has been divided into mutually-exclusive compartments comprising of individuals that are susceptible ( $S(a,t)$ ), individuals exposed to HSV-2 but show no clinical symptoms of the disease ( $E(a,t)$ ), individuals with infection (virus shedding) with clinical symptoms of HSV-2 ( $H(a,t)$ ), and individuals with infection, which is quiescent (an infection is present but not active in a host) ( $Q(a,t)$ ).

$\displaystyle N(a,t)=S(a,t)+E(a,t)+H(a,t)+Q(a,t).$

Eventually, the population density evolves according to the following equation:

$\displaystyle\frac{\partial N(a,t)}{\partial a}+\frac{\partial N(a,t)}{% \partial t}=-\mu(a)N(a,t),$

where $\mu(a)$ represents the instantaneous age-dependent death rate (the probability that an individual of age $a$ dies before reaching the age $a+\Delta a$ is $\mu(a)\Delta a+O(\Delta a))$ .

Following [21], it is assumed that the net reproductive rate of the population equals to unity, and the total population is at an equilibrium.

$\displaystyle\int_{0}^{A}v(a)e^{-\int_{0}^{a}\mu(\tau)d\tau}da=1;\quad N(a,t)=% N_{\infty}(a)=v_{0}e^{-\int_{0}^{a}\mu(\tau)d\tau},$ (1)

where $v(a)$ is the age-dependent fertility rate, $N_{\infty}(a)$ is the steady-state population distribution and $v_{0}$ is the steady-state number of new-borns.

Let $\beta(a,b)$ be the age-dependent contact rate. In other words, the probability of contact between a susceptible of age $a$ and an infectious individual of age $b$ during the time interval $[t,t+\Delta t]$ is given by $\beta(a,b)\Delta t+O(\Delta t)$ . Hence, the associated force of infection (i.e., the rate at which new cases of infection are generated) at time $t$ , denoted by $\lambda(a,t)$ , is given by (this form of the force of infection has been used in [13, 20]):

$\displaystyle\lambda(a,t)=\int_{0}^{A}\beta(a,b)[H(b,t)+Q(b,t)]db.$

Figure 1.

Flow diagram of the model Eq. (2).

The age-structured model for the transmission dynamics of HSV-2 in a population is given by the following system of partial differential equations:

$\displaystyle\frac{\partial S}{\partial a}+\frac{\partial S}{\partial t}=-[\mu% (a)+\lambda(a,t)]S(a,t),$ $\displaystyle\frac{\partial E}{\partial a}+\frac{\partial E}{\partial t}=% \lambda(a,t)S(a,t)-[\mu(a)+\sigma]E(a,t),$ (2) $\displaystyle\frac{\partial H}{\partial a}+\frac{\partial H}{\partial t}=% \sigma E(a,t)+rQ(a,t)-[\mu(a)+q]H(a,t),$ $\displaystyle\frac{\partial Q}{\partial a}+\frac{\partial Q}{\partial t}=qH(a,% t)-rQ(a,t)-\mu(a)Q(a,t),$

subject to the boundary and initial conditions:

$\displaystyle S(0,t)=v_{0},\quad E(0,t)=0,\quad H(0,t)=0,\quad Q(0,t)=0,$ (3) $\displaystyle S(a,0)=S_{0}(a),\quad E(a,0)=E_{0}(a),\quad H(a,0)=H_{0}(a),% \quad Q(a,0)=Q_{0}(a).$

In Eq. (2), $\mu(a)$ is the age-dependent natural death rate, $\sigma$ is the rate at which latent individuals develop symptoms, $r$ is the re-activation rate of symptoms by individuals in the quiescent state, $q$ is the rate at which infectious individuals revert to quiescent state. Unlike [23], it is assumed that disease-induced mortality is negligible (so that infected individuals do not die of the disease) and that the quiescent infectious individuals can transmit infection at the same rate as the non-quiescent infectious individuals. A flow diagram of the model is depicted in Fig. 1, and the associated parameters are described in Table 1.

The model extends the model in [23] by considering the age-densities of all individuals.

Furthermore, the model Eq. (2) is an extension of the HSV-2 models in [2, 6], by including classes for latent ( $E$ ).

2.1 Basic properties

Consider the following change of variables

$\displaystyle x_{1}(a,t)=\frac{S(a,t)}{N_{\infty}(a)},\quad x_{2}(a,t)=\frac{E% (a,t)}{N_{\infty}(a)},\quad x_{3}(a,t)=\frac{H(a,t)+Q(a,t)}{N_{\infty}(a)}.$ (4)

Table 1

Description of variables and parameters of the model Eq. (2)

Variable	Description
$S(t)$	Population of susceptible individuals
$E(t)$	Population of exposed individuals
$H(t)$	Population of infectious (symptomatic) individuals
$Q(t)$	Population of quiescent infectious individuals
Parameter	Description
$\mu(a)$	Age-dependent natural death rate
$\beta(a,b)$	Age-dependent Effective contact rate
$\sigma$	Progression rate from exposed to infectious class
$q$	Rate at which infectious individuals revert to quiescent state
$r$	Activation rate of infectious individuals in the quiescent state

Adding the last two equations of Eq. (2) and then using Eq. (4) in Eq. (2) gives the following transformed system [25]:

$\displaystyle\frac{\partial x_{1}}{\partial a}+\frac{\partial x_{1}}{\partial t% }=-\lambda_{1}(a,t)x_{1}(a,t),$ $\displaystyle\frac{\partial x_{2}}{\partial a}+\frac{\partial x_{2}}{\partial t% }=\lambda_{1}(a,t)x_{1}(a,t)-\sigma x_{2}(a,t),$ (5) $\displaystyle\frac{\partial x_{3}}{\partial a}+\frac{\partial x_{3}}{\partial t% }=\sigma x_{2}(a,t),$

and associated boundary and initial conditions:

$\displaystyle x_{1}(0,t)=1,\quad x_{2}(0,t)=0,\quad x_{3}(0,t)=0,$ (6) $\displaystyle x_{1}(a,0)=f_{1}(a),\quad x_{2}(a,0)=f_{2}(a),\quad x_{3}(a,0)=f% _{3}(a).$

The force of infection, $\lambda(a,t)$ , now becomes

$\displaystyle\lambda_{1}(a,t)=\int_{0}^{A}\beta(a,b)x_{3}(b,t)N_{\infty}(b)db.$ (7)

Now, introducing a new variable $\bar{x}_{1}(a,t)=x_{1}(a,t)-1$ into the transformed system Eq. (2.1) gives

$\displaystyle\frac{\partial\bar{x}_{1}}{\partial a}+\frac{\partial\bar{x}_{1}}% {\partial t}=-\lambda_{1}(a,t)\bar{x}_{1}(a,t),$ $\displaystyle\frac{\partial x_{2}}{\partial a}+\frac{\partial x_{2}}{\partial t% }=\lambda_{1}(a,t)\bar{x}_{1}(a,t)-\sigma x_{2}(a,t),$ (8) $\displaystyle\frac{\partial x_{3}}{\partial a}+\frac{\partial x_{3}}{\partial t% }=\sigma x_{2}(a,t),$

and associated boundary and initial conditions:

$\displaystyle\bar{x}_{1}(0,t)=0,\quad x_{2}(0,t)=0,\quad x_{3}(0,t)=0,$ (9) $\displaystyle\bar{x}_{1}(a,0)=f_{1}(a)-1,\quad x_{2}(a,0)=f_{2}(a),\quad x_{3}% (a,0)=f_{3}(a).$

Re-write the initial-boundary-value problem Eqs (2.1) and (2.1) as an abstract Cauchy problem as follows [13]: Let the Banach space,

$\displaystyle X=L^{1}(0,A)\times L^{1}(0,A)\times L^{1}(0,A)$

with the norm $\|\phi\|=\sum_{j=1}^{3}\|\phi_{j}\|_{1}$ , where $\phi(a)=(\phi_{1}(a),\phi_{2}(a),\phi_{3}(a))^{T}\in X$ and $\|\phi_{j}\|_{1}$ is the norm of $L^{1}(0,A)$ ) by defining a linear operator, $\mathfrak{A}$ , as follows:

$\displaystyle(\mathfrak{A}\phi)(a)=\left(\begin{array}[]{c}-\frac{d\phi_{1}(a)% }{da}\\ -\frac{d\phi_{2}(a)}{da}-\sigma\phi_{2}\\ 0\end{array}\right),$ (10)

with $\phi(a)\in D(\mathfrak{A})=\{\phi\in X:\phi_{j}\text{ is absolutely continuous% on }[0,A],\phi(0)=0\}$ .

Let $F:X\rightarrow X$ be the non-linear operator given by (consisting of the remaining terms, excluding those used in the definition of the operator $A$ , in Eq. (2.1))

$\displaystyle(F\phi)(a)=\left(\begin{array}[]{c}-(P(\phi_{3}))(a)(1+\phi_{1})% \\ (P(\phi_{3}))(a)(1+\phi_{1})\\ \sigma\phi_{2}(a)\end{array}\right),\quad\phi\in X$ (11)

where $P:L^{1}(0,A)\rightarrow L^{\infty}(0,A)$ is a linear bounded operator, given by

$\displaystyle(P(h)(a)=\int_{0}^{A}\beta(a,b)N_{\infty}(b)h(b)db.$

Thus, it follows from the above definitions that the initial-boundary-value problem Eqs (2.1) and (2.1) can be re-written as the following abstract Cauchy problem in $X$ :

$\displaystyle\frac{du(t)}{dt}=\mathfrak{A}u(t)+F(u(t)),\quad u(0)=u_{0}\in X,$ (12)

where $u(t)=(\bar{x}_{1}(\cdot,t),x_{2}(\cdot,t),x_{3}(\cdot,t))^{T}$ and $u_{0}=(f_{1}(a)-1,f_{2}(a),f_{3}(a))^{T}$ . Furthermore, it follows from Eqs (10) and (11) that the operator $\mathfrak{A}$ is the infinitesimal generator of $C_{0}$ -semigroup $T(t),t\geqslant$ 0 and $F$ is continuously Frechet differentiable on $X$ [20]. Thus, using Proposition 4.16 of [39], the results below are established.

Theorem 1 (i) For each $u_{0}\in X$ , there exists a maximal interval of existence $[0,b)$ and unique continuous mild solution $t\rightarrow u(t,u_{0})$ such that:

$\displaystyle u(t,u_{0})=T(t)u_{0}+\int_{0}^{t}T(t-s)F(u(s,u_{0}))ds\quad\text% { for all }t\in[0,b),$

with either $b=\infty$ or $\lim\limits_{t\rightarrow b^{-}}\|u(t,u_{0})\|=\infty$ .

(ii) If $u_{0}\in D(\mathfrak{A})$ , then $u(t,u_{0})\in D(\mathfrak{A})$ for $t\in[0,b)$ and the mild solution $u(t,u_{0})$ is continuously differentiable and satisfies Eq. (12) on the interval $[0,b)$ .

Define,

$\displaystyle\mathcal{D}=\{(\bar{x}_{1},x_{2},x_{3})\in X:\bar{x}_{1}\geqslant% -1,x_{2}\geqslant 0,x_{3}\geqslant 0\},$

and,

$\displaystyle\mathcal{D}_{0}=\{(\bar{x}_{1},x_{2},x_{3})\in X:-1\leqslant\bar{% x}_{1}\leqslant 0,0\leqslant x_{2}\leqslant 1,0\leqslant x_{3}\leqslant 1\}.$

The following result is established.

Lemma 1 The mild solution of the abstract Cauchy problem Eq. (12), given by $u(t,u_{0})$ with $u_{0}\in\mathcal{D}$ , enters $\mathcal{D}_{0}$ in finite time, and the set $\mathcal{D}_{0}$ is positively-invariant.

Proof..

Following [13, 20]. First of all, the solution of the system Eq. (2.1), along the characteristic lines, can be represented [10, 39, 31] by fixing arbitrary time $(t_{1})$ and age $(a_{1})$ and considering the functions $\tilde{x}_{1}(h)=x_{1}(a_{1}+h,t_{1}+h)$ and $\tilde{\lambda}(h)=\lambda_{1}(a_{1}+h,t_{1}+h)$ . Thus, the first equation of the system Eq. (2.1), becomes

$\displaystyle\frac{d\tilde{x}_{1}}{dh}+\tilde{\lambda}\tilde{x}_{1}=0,$ (13)

The solution of Eq. (13)

$\displaystyle\tilde{x}_{1}(h)=\tilde{x}_{1}(0)\exp\left(-\int_{0}^{h}\tilde{% \lambda}(\tau)d\tau\right).$ (14)

Thus,

$\displaystyle x_{1}(a_{1}+h,t_{1}+h)=x_{1}(a_{1},t_{1})\exp\left(-\int_{0}^{h}% \lambda_{1}(a_{1}+\tau,t_{1}+\tau)d\tau\right).$

Setting $(a_{1},t_{1})=(a-t,0)$ and $h=t$ for $a>t$ in Eq. (14), it follows that

$\displaystyle x_{1}(a,t)=f_{1}(a-t)\exp\left(-\int_{0}^{t}\lambda_{1}(a-t+\tau% ,\tau)d\tau\right).$

Similarly, for the case when $a<t$ , setting $(a_{1},t_{1})=(0,t-a)$ and $h=a$ in Eq. (14) gives

$\displaystyle x_{1}(a,t)=\exp\left(-\int_{0}^{a}\lambda_{1}(\tau,t-a+\tau)d% \tau\right).$

Hence,

$\displaystyle x_{1}(a,t)=\left\{\begin{array}[]{ll}f_{1}(a-t)\exp\left(-\int_{% 0}^{t}\lambda_{1}(a-t+\tau,\tau)d\tau\right),&a>t\\ \exp\left(-\int_{0}^{a}\lambda_{1}(\tau,t-a+\tau)d\tau\right),&t>a\end{array}\right.$ (15)

It is clear that $x_{1}(a,t)\geqslant$ 0 whenever $f_{1}(a)\geqslant$ 0. In similar way the last equation of Eq. (2.1) gives

$\displaystyle x_{3}(a,t)=\left\{\begin{array}[]{ll}f_{3}(a-t)+\sigma\int_{0}^{% t}x_{2}(a-t+\tau,\tau)d\tau,&a>t\\ \sigma\int_{0}^{a}x_{2}(\tau,t-a+\tau)d\tau,&t>a,\end{array}\right.$ (16)

Substituting Eq. (16) into Eq. (7) gives

$\displaystyle\lambda_{1}(t)=(Gx_{2})(a,t),$ (17)

where $G$ is the transformation from $x_{2}(a,t)$ to $\lambda_{1}(t)$ . The second equation of the system Eq. (2.1) can be written as an abstract Cauchy problem:

$\displaystyle\frac{dx_{2}}{dt}=Bx_{2}(a,t)+(Gx_{2}(a,t))(x_{1}),\quad x_{2}(0)% =f_{2},$

where $B$ is an operator given by

$\displaystyle B=-\frac{d}{da}-(\kappa+\sigma),\text{ with }D(B)=D(\mathfrak{A}).$

Hence, $x_{2}$ can be found as follows:

$\displaystyle x_{2}(t)=M(t)f_{2}(0)+\int_{0}^{t}M(t-\tau)(Gx_{2}(\tau))% \mathfrak{(\tau)}d\tau,$ (18)

where, $M(t)=\exp{(tB)}$ is the positive $C_{0}$ -semigroup generated by the closed operator $B$ [13]. If $f_{2}\geqslant$ 0 and $f_{3}\geqslant$ 0, then $G$ and $M(t),t\geqslant$ 0 are positive. Thus, Eq. (18) shows that $x_{2}$ is positive. Since $x_{2}$ is positive, it follows that $x_{3}$ is positive. Hence, $u(t,u_{0})\in\mathcal{D}$ for all $t\geqslant$ 0 whenever $u_{0}\in\mathcal{D}_{0}$ . ∎

It follows from Lemma 1 that the norm of the local solution, $u(t,u_{0}),u_{0}\in D(\mathfrak{A})\cap\mathcal{D}$ , of the abstract Cauchy problem Eq. (12) is finite as long as it is defined. Hence, the following result is established:

Theorem 2 The abstract Cauchy problem Eq. (12) has a unique global classical solution on $X$ with respect to initial data $u_{0}\in D(\mathfrak{A})\cap\mathcal{D}$ .

Theorem 2 shows that the initial-boundary-value problem Eqs (2.1) and (2.1) (or, equivalently, Eqs (2) and (2)) has a unique positive global solution with respect to the positive initial data.

3. Disease-free equilibrium (DFE)

Consider the following separable contact rate [21, 22, 12]:

$\displaystyle\beta(a,b)=\beta_{1}(a)\beta_{2}(b),$ (19)

then the analysis of the initial-boundary-value problem Eqs (2.1) and (2.1) will be explored when the contact rate is separable. Using the separable contact rate Eq. (19) in Eq. (7) gives:

$\displaystyle\lambda_{2}(a,t)=\beta_{1}(a)\int_{0}^{A}\beta_{2}(b)x_{3}(b,t)N_% {\infty}(b)db.$ (20)

3.1 Local stability of DFE

The DFE of the system Eq. (2.1) is given by

$\displaystyle\mathcal{E}_{0}=(x_{1},x_{2},x_{3})=(1,0,0).$

Following [21, 22], then the local stability of the DFE is studied by considering exponential solutions of the form:

$\displaystyle x_{1}(a,t)=1+\bar{x}_{1}(a)e^{\zeta t},\quad x_{2}(a,t)=\bar{x}_% {2}(a)e^{\zeta t},\quad x_{3}(a,t)=\bar{x}_{3}(a)e^{\zeta t},$ (21)

Linearizing Eq. (2.1) about the DFE ( $\mathcal{E}_{0}$ ) gives,

$\displaystyle\zeta\bar{x}_{2}(a)+\frac{\bar{x}_{2}(a)}{da}=\lambda_{2}(a)-% \sigma\bar{x}_{2}(a),$ (22) $\displaystyle\zeta\bar{x}_{3}(a)+\frac{\bar{x}_{3}(a)}{da}=\sigma\bar{x}_{2}(a),$

where,

$\displaystyle\lambda_{2}(a)=\beta_{1}(a)\int_{0}^{A}\beta_{2}(b)\bar{x}_{3}(b)% N_{\infty}(b)db=\beta_{1}(a)V_{0},$

with

$\displaystyle V_{0}=\int_{0}^{A}\beta_{2}(b)\bar{x}_{3}(b)N_{\infty}(b)db;% \quad V_{0}\neq 0.$ (23)

Solving system Eq. (22) gives,

$\displaystyle\bar{x}_{2}(a)=V_{0}\int_{0}^{a}\beta_{1}(u)e^{-(\zeta+\sigma)(a-% u)}du,$ (24) $\displaystyle\bar{x}_{3}(a)=V_{0}\sigma\int_{0}^{a}\int_{0}^{u}\beta_{1}(u_{1}% )e^{-(\zeta+\sigma)(u-u_{1})}e^{-\zeta(a-u)}du_{1}du=V_{0}Q_{1}(a,\zeta),$

where,

$\displaystyle Q_{1}(a,\zeta)=\sigma\int_{0}^{a}\int_{0}^{u}\beta_{1}(u_{1})e^{% -(\zeta+\sigma)(u-u_{1})}e^{-\zeta(a-u)}du_{1}du.$

Substituting Eq. (24) into Eq. (23) gives,

$\displaystyle V_{0}=V_{0}\int_{0}^{A}\beta_{2}(b)Q_{1}(b,\zeta)N_{\infty}(b)db,$ (25)

so that (by dividing both sides of the equation Eq. (25) by $V_{0}$ )

$\displaystyle 1=\int_{0}^{A}\beta_{2}(b)Q_{1}(b,\zeta)N_{\infty}(b)db.$ (26)

Let,

$\displaystyle\mathcal{G}(\zeta)=\int_{0}^{A}\beta_{2}(b)Q_{1}(b,\zeta)N_{% \infty}(b)db.$ (27)

Define the reproduction number as follows:

$\displaystyle\mathcal{R}_{0}=\int_{0}^{A}\beta_{2}(b)Q_{1}(b,0)N_{\infty}(b)db.$

The function $\mathcal{G}$ in Eq. (27) satisfies the following properties:

(i)

$\mathcal{G}(0)=\mathcal{R}_{0}$ ;

(ii)

$\mathcal{G}$ is a monotone decreasing function in $\zeta$ ;

(iii)

$\lim\limits_{\zeta\rightarrow\infty}\mathcal{G}(\zeta)=0$ .

Using the above properties of the function $\mathcal{G}$ and Eq. (26) the following result is established.

Theorem 3 The DFE of the model Eq. (2.1) with Eq. (20), given by $\mathcal{E}_{0}$ , is locally-asymptotically stable (LAS) if $\mathcal{R}_{0}<$ 1, and unstable if $\mathcal{R}_{0}>$ 1.

Theorem 3 shows that the spread of HSV-2 can be effectively controlled if the initial sizes of the sub-populations of the model are in the basin of attraction of the DFE ( $\mathcal{E}_{0}$ ). For effective disease elimination to be independent of such initial sizes, a global asymptotic stability result has to be proved for the DFE. This is explored below.

3.2 Global stability of DFE

We claim the following result.

Theorem 4 The DFE of the model Eq. (2.1) with Eq. (20) is globally-asymptotically stable (GAS) in $\mathcal{D}_{0}$ if $\mathcal{R}_{0}<$ 1.

Proof..

The proof is similar to the proof of Theorem 4 given in [25]. Following [21, 22], let

$\displaystyle f(a,t)=\lambda_{2}(a,t)x_{1}(a,t).$ (28)

Noting that $x_{1}(a,t)\leqslant$ 1 in $\mathcal{D}_{0}$ , it follows that

$\displaystyle f(a,t)\leqslant\lambda_{2}(a,t).$ (29)

Integrating Eq. (2.1), with Eq. (20), along the characteristic lines whenever $a<t$ , gives,

$\displaystyle x_{2}(a,t)=\int_{0}^{a}e^{-\sigma(a-u)}f(u,t)du,$ (30) $\displaystyle x_{3}(a,t)=\sigma\int_{0}^{a}x_{2}(u,t)du=\sigma\int_{0}^{a}\int% _{0}^{u}e^{-\sigma(u-u_{1})}f(u_{1},t)du_{1}du,$

Substituting Eq. (30) in Eq. (29) gives,

$\displaystyle f(a,t)\leqslant\beta_{1}(a)V(t)=\beta_{1}(a)\int_{0}^{A}N_{% \infty}(b)\beta_{2}(b)\bigg{[}\sigma\int_{0}^{a}\int_{0}^{u}e^{-\sigma(u-u_{1}% )}f(u_{1},t)du_{1}du\bigg{]}db.$ (31)

Let,

$\displaystyle F(a)=\limsup_{t\rightarrow\infty}{f(a,t)}.$

Taking the $\limsup$ (when $t\rightarrow\infty$ ) of both sides of Eq. (31), and using Fatou’s Lemma [24] gives,

$\displaystyle F(a)\leqslant\beta_{1}(a)L,$ (32)

where,

$\displaystyle L=\int_{0}^{A}N_{\infty}(b)\beta_{2}(b)\bigg{[}\sigma\int_{0}^{a% }\int_{0}^{u}e^{-\sigma(u-u_{1})}F(u_{1})du_{1}du\bigg{]}db.$ (33)

Using inequality Eq. (32) in Eq. (34) gives

$\displaystyle L\leqslant\int_{0}^{A}N_{\infty}(b)\beta_{2}(b)\bigg{[}\sigma% \int_{0}^{a}\int_{0}^{u}e^{-\sigma(u-u_{1})}\beta(u_{1})Ldu_{1}du\bigg{]}db=L% \mathcal{R}_{0}.$ (34)

It follows that whenever $\mathcal{R}_{0}<$ 1, then $L=$ 0. Thus,

$\displaystyle\limsup_{t\rightarrow\infty}{f(a,t)}=F(a)=0.$

Hence (from Eq. (28)),

$\displaystyle\limsup_{t\rightarrow\infty}{\lambda_{2}(a,t)}=0.$ (35)

It follows then (by using Eq. (35) in Eq. (20)) that

$\displaystyle\lim_{t\rightarrow\infty}x_{3}(a,t)=0.$

By using comparison theorem [29] it can be shown that

$\displaystyle\lim_{t\rightarrow\infty}x_{2}(a,t)=0.$

Finally, using the relation $x_{1}(a,t)=1-x_{2}(a,t)-x_{3}(a,t)$ , gives

$\displaystyle\lim_{t\rightarrow\infty}x_{1}(a,t)=1.$

Hence, $\lim_{t\rightarrow\infty}(x_{1}(a,t),x_{2}(a,t),x_{3}(a,t))=(1,0,0)$ , as required. ∎

4. Endemic equilibria

4.1 Existence and uniqueness

The existence of endemic equilibria of the initial-boundary-value problem Eqs (2.1) and (2.1) (which is equivalent to Eqs (2.1) and (2.1)) will be explored. Let,

$\displaystyle\mathcal{E}=(x_{1}^{*}(a),x_{2}^{*}(a),x_{3}^{*}(a))$

represents any arbitrary equilibrium point of the model Eqs (2.1) and (2.1) with Eq. (20). The associated force of infection Eq. (20) at steady-state as is given by

$\displaystyle\lambda_{2}^{*}(a)=\beta_{1}(a)\int_{0}^{A}\beta_{2}(b)x_{3}^{*}(% b)N_{\infty}(b)db=\beta_{1}(a)V^{*},$ (36)

where,

$\displaystyle V^{*}=\int_{0}^{A}\beta_{2}(b)x_{3}^{*}(b)N_{\infty}(b)db.$ (37)

Solving the equations of the system Eq. (2.1) at steady-state, and noting Eq. (36), gives

$\displaystyle x_{1}^{*}(a)=e^{-V^{*}\int_{0}^{a}\beta_{1}(\tau)d\tau},$ $\displaystyle x_{2}^{*}(a)=V^{*}e^{-\sigma a}\int_{0}^{a}\beta_{1}(u_{1})e^{-V% ^{*}\int_{0}^{u_{1}}\beta_{1}(\tau)d\tau}e^{\sigma u_{1}}du_{1},$ (38) $\displaystyle x_{3}^{*}(a)=V^{*}\sigma\int_{0}^{a}e^{-\sigma u_{1}}\int_{0}^{u% _{1}}\beta_{1}(u_{2})e^{-V^{*}\int_{0}^{u_{2}}\beta_{1}(\tau)d\tau}e^{\sigma u% _{2}}du_{2}du_{1}=V^{*}M_{1}(a,V^{*}),$

where,

$\displaystyle M_{1}(a,V^{*})=\sigma\int_{0}^{a}e^{-\sigma u_{1}}\int_{0}^{u_{1% }}\beta_{1}(u_{2})e^{-V^{*}\int_{0}^{u_{2}}\beta_{1}(\tau)d\tau}e^{\sigma u_{2% }}du_{2}du_{1}.$

Substituting the expressions for Eq. (4.1) in Eq. (37) gives,

$\displaystyle V^{*}=V^{*}\int_{0}^{A}\beta_{2}(a)M_{1}(a,V^{*})N_{\infty}(a)da.$ (39)

It follows from equation Eq. (39) that the quantity $V^{*}$ is either 0 or satisfies the following equation:

$\displaystyle 1=\int_{0}^{A}\beta_{2}(a)M_{1}(a,V^{*})N_{\infty}(a)da.$ (40)

Let,

$\displaystyle\mathcal{H}(V^{*})=\int_{0}^{A}\beta_{2}(a)M_{1}(a,V^{*})N_{% \infty}(a)da.$

Since $\mathcal{H}$ satisfies

(i)

$\mathcal{H}(0)=\mathcal{R}_{0}$ ;

(ii)

$\mathcal{H}$ is a monotone decreasing function in $V^{*}$ ;

(iii)

$\lim\limits_{V^{*}\rightarrow\infty}\mathcal{H}(V^{*})=0$ ,

then, the result below can be established.

Theorem 5 The model Eq. (2.1) with Eq. (20) has a unique endemic equilibrium point (EEP), given by Eq. (4.1), whenever $\mathcal{R}_{0}>$ 1.

4.2 Local stability of EEP

Following [25], the local stability of the unique endemic equilibrium Eq. (4.1) is now explored as follows. Let $\hat{x}_{1},\hat{x}_{2},\hat{x}_{3}$ be the perturbations of ${x}_{1}^{*},{x}_{2}^{*},{x}_{3}^{*}$ .

Consider, the exponential solutions:

$\displaystyle\hat{x}_{1}(a,t)=\bar{x}_{1}(a)e^{\omega t},\quad\hat{x}_{2}(a,t)% =\bar{x}_{2}(a)e^{\omega t},\quad\hat{x}_{3}(a,t)=\bar{x}_{3}(a)e^{\omega t}$ (41)

Using Eq. (41) in Eq. (2.1) with Eq. (36) gives the following system:

$\displaystyle\frac{d\bar{x}_{1}(a)}{da}=-\omega\bar{x}_{1}(a)-\beta_{1}(a)V^{*% }(a)\bar{x}_{1}(a)-\beta_{1}(a)x_{1}^{*}(a)\bar{V}(a)$ $\displaystyle\frac{d\bar{x}_{2}(a)}{da}=-\omega\bar{x}_{2}(a)+\beta_{1}(a)V^{*% }(a)\bar{S}(a)+\beta_{1}(a)x_{1}^{*}(a)\bar{V}(a)-\sigma\bar{x}_{2}(a)$ (42) $\displaystyle\frac{d\bar{x}_{3}(a)}{da}=-\omega\bar{x}_{3}(a)+\sigma\bar{x}_{2% }(a)$

with initial conditions $\hat{x}_{1}(0)=0,\hat{x}_{2}(0)=0,\hat{x}_{3}(0)=0$ , where,

$\displaystyle\bar{V}(a)=\int_{0}^{A}\beta_{2}(b)\bar{x}_{3}N_{\infty}(b)db.$ (43)

Consider the following new variables:

$\displaystyle y_{1}(a)=\frac{\bar{x}_{1}(a)}{\bar{V}(a)},y_{2}(a)=\frac{\bar{x% }_{2}(a)}{\bar{V}(a)},y_{3}(a)=\frac{\bar{x}_{3}(a)}{\bar{V}(a)}.$ (44)

Using Eq. (44) into the system Eq. (4.2) gives,

$\displaystyle\frac{dy_{1}(a)}{da}=-\omega y_{1}-\beta_{1}(a)V^{*}y_{1}-\beta_{% 1}(a)x_{1}^{*},$ $\displaystyle\frac{dy_{2}(a)}{da}=-\omega y_{2}+\beta_{1}(a)V^{*}y_{1}+\beta_{% 1}(a)x_{1}^{*}-\sigma y_{2},$ (45) $\displaystyle\frac{dy_{3}(a)}{da}=-\omega y_{3}+\sigma y_{2},$

$\displaystyle 1=\int_{0}^{A}\beta_{2}(a)y_{3}(a)N_{\infty}(a)da.$ (46)

Solving the above system gives

$\displaystyle y_{1}(a)=\int_{0}^{a}{-\beta_{1}(u)x_{1}^{*}(u)e^{-\omega(a-u)}e% ^{\int_{u}^{a}{\beta_{1}(u_{1})V^{*}(u_{1})}du_{1}}}du,$ $\displaystyle y_{2}(a)=\int_{0}^{a}[{\beta_{1}(u)x_{1}^{*}(u)+\beta_{1}(u)V^{*% }(u)y_{1}(u)}]e^{-(\omega+\sigma)(a-u)}du,$ (47) $\displaystyle y_{3}(a)=\int_{0}^{a}\sigma y_{2}(u)e^{-(\omega)(a-u)}du.$

Substituting the equations in Eq. (4.2) into Eq. (46) gives

$\displaystyle 1=\int_{0}^{A}\beta_{2}(a)y_{3}(a)N_{\infty}(a)da.$ (48)

Define,

$\displaystyle\mathcal{Q}(\omega)=\int_{0}^{A}\beta_{2}(a)y_{3}(a)N_{\infty}(a)da.$ (49)

It’s clear that the function $\mathcal{Q}(\omega)$ is a monotone decreasing function in $\omega$ and $\lim\limits_{\omega\rightarrow\infty}\mathcal{Q}(\omega)=0$ .

Now consider $\mathcal{P}(\omega)=\mathcal{Q}(\omega)-1$ , it is clear that if $\mathcal{Q}(0)<$ 1, then $\mathcal{P}(\omega)$ has a negative root. Thus the following result is established.

Theorem 6 The unique endemic equilibrium of the model Eq. (2.1) with Eq. (20) is LAS whenever $\mathcal{R}_{0}>$ 1 and $\mathcal{Q}(0)<$ 1.

The consequence of Theorem 6 is that the disease will persist (become endemic) in the community whenever $\mathcal{R}_{0}>1$ and $\mathcal{Q}(0)<1$ . In summary, it is shown that the model Eq. (2.1) with Eq. (20) has a globally-asymptotically stable DFE in $\mathcal{D}_{0}$ whenever the associated reproduction number $\mathcal{R}_{0}$ is less than unity. Furthermore, the model has a unique endemic equilibrium whenever the reproduction number exceeds unity. This equilibrium is shown to be locally-asymptotically stable if another condition holds ( $\mathcal{Q}(0)<$ 1). The global stability of the EEP is now considered for another special case, as described below.

5. Model with constant natural death and contact rates

Consider a special case of the model Eq. (2) where the natural death rate is constant (i.e., $\mu$ is independent of age) and the contact rate ( $\beta(a,b)$ ) is constant (denoted by $\beta(a,b)=\beta$ ; see also [7, 8, 15]). Consider, now, the model Eq. (2) with $\mu(a)$ replaced by the constant $\mu$ and $\beta(a,b)$ by a constant $\beta$ (so that the force of infection, $\lambda(a,t)$ , is now given by $\lambda(t)=\int_{0}^{A}\beta H(a,t)+Q(a,t)da$ ), given by:

$\displaystyle\frac{\partial S}{\partial a}+\frac{\partial S}{\partial t}=-[\mu% (a)+\lambda(a,t)]S(a,t),$ $\displaystyle\frac{\partial E}{\partial a}+\frac{\partial E}{\partial t}=% \lambda(a,t)S(a,t)-[\mu(a)+\sigma]E(a,t),$ (50) $\displaystyle\frac{\partial H}{\partial a}+\frac{\partial H}{\partial t}=% \sigma E(a,t)+rQ(a,t)-[\mu(a)+q]H(a,t),$ $\displaystyle\frac{\partial Q}{\partial a}+\frac{\partial Q}{\partial t}=qH(a,% t)-rQ(a,t)-\mu(a)Q(a,t),$

subject to the boundary and initial conditions:

$\displaystyle S(0,t)=\mu N_{0},\quad E(0,t)=0,\quad H(0,t)=0,\quad Q(0,t)=0,$ (51) $\displaystyle S(a,0)=S_{0}(a),\quad E(a,0)=E_{0}(a),\quad H(a,0)=H_{0}(a),% \quad Q(a,0)=Q_{0}(a).$

Let,

$\displaystyle\widetilde{S}(t)=\int_{0}^{A}S(a,t)da,\widetilde{E}(t)=\int_{0}^{% A}E(a,t)da,$ $\displaystyle\widetilde{H}(t)=\int_{0}^{H}Q(a,t)da,\widetilde{Q}(t)=\int_{0}^{% A}Q(a,t)da,$

and consider the following assumptions in line with [7]:

(A1)

$\widetilde{S}(t),\widetilde{E}(t),\widetilde{H}(t),\widetilde{Y}(t)$ , and $\widetilde{Q}(t)$ are finite with $\widetilde{S}(t)+\widetilde{E}(t)+\widetilde{H}(t)+\widetilde{Q}(t)=N_{0}$ ;

(A2)

The equation

$\displaystyle\frac{d}{dt}X(t)=\int_{0}^{A}{\frac{\partial X}{\partial t}}da$

is valid for all $X=\widetilde{S},\widetilde{E},\widetilde{H}$ , and $\widetilde{Q}$ ;

(A3)

Each of $S(a,0),E(a,0),H(a,0)$ , and $Q(a,0)$ is positive, continuous, differentiable and tends to zero as $a$ tends to $A$ ;

(A4)

$\lim\limits_{a\rightarrow 0}X(a,0)=\lim_{t\rightarrow 0}X(0,t)$ and the right derivative of $X(a,0)=$ 0 at $a=$ 0, for all $X={S},{E},{H}$ , and ${Q}$ .

By Assumptions A1–A4, integrating the set of partial differential equations in Eq. (5) with respect to $a$ , gives the following system of ordinary differential equations:

$\displaystyle\frac{d\widetilde{S}}{dt}=\mu N_{0}-[\mu+\widetilde{\lambda}(t)]% \widetilde{S}(t),$ $\displaystyle\frac{d\widetilde{E}}{dt}=\widetilde{\lambda}(t)\widetilde{S}(t)-% (\mu+\sigma)\widetilde{E}(t),$ (52) $\displaystyle\frac{d\widetilde{H}}{dt}=\sigma\widetilde{E}(t)+r\widetilde{Q}-(% q+\mu)\widetilde{H}(t),$ $\displaystyle\frac{d\widetilde{Q}}{dt}=q\widetilde{H}(t)-(\mu+r)\widetilde{Q}(% t),$

adding the last equations of Eq. (5) gives

$\displaystyle\frac{d\widetilde{S}}{dt}=\mu N_{0}-[\mu+\widetilde{\lambda}(t)]% \widetilde{S}(t),$ $\displaystyle\frac{d\widetilde{E}}{dt}=\widetilde{\lambda}(t)\widetilde{S}(t)-% (\mu+\sigma)\widetilde{E}(t),$ (53) $\displaystyle\frac{d\widetilde{I}}{dt}=\sigma\widetilde{E}(t)-\mu\widetilde{I}% (t),$

$\displaystyle\widetilde{\lambda}(t)=\beta\widetilde{I}.$ (54)

Define the following positively-invariant region for the model Eq. (5):

$\displaystyle\mathcal{D}_{c}=\biggl{\{}(\widetilde{S},\widetilde{E},\widetilde% {I})\in\mathbb{R}^{3}_{+}:\widetilde{S}+\widetilde{E}+\widetilde{I}\leqslant N% _{0}\biggr{\}}.$

and the associated reproduction number of the model Eq. (5) is given by

$\displaystyle\mathcal{R}_{cons}=\frac{\beta N_{0}}{\mu},$ (55)

Let $\mathcal{E}_{1c}=(S^{**},E^{**},I^{**})$ be an endemic equilibrium point of the model Eq. (5). It should be mentioned that since the model Eq. (5) is a special case of the model Eq. (2), then the following results are established

Corollary 1 The DFE of the model Eq. (5) is globally-asymptotically stable whenever $\mathcal{R}_{cons}<$ 1.

Corollary 2 The model Eq. (5) has a unique endemic equilibrium point, denoted by $\mathcal{E}_{1c}$ , whenever $\mathcal{R}_{cons}>$ 1.

It is convenient to define the region

$\displaystyle\mathcal{D}_{01}=\biggl{\{}(\widetilde{S},\widetilde{E},% \widetilde{I})\in\mathcal{D}_{c}:\widetilde{S}=N_{0},\widetilde{E}=\widetilde{% I}=0\biggr{\}}.$

The following result can be established.

Theorem 7 The unique endemic equilibrium of the model Eq. (5) is GAS in $\mathcal{D}_{c}\setminus\mathcal{D}_{01}$ if $\mathcal{R}_{cons}>$ 1.

Proof..

Consider the following Lyapunov function

$\displaystyle\mathcal{F}=\frac{1}{2}\left[(\widetilde{S}-S^{**})+(\widetilde{E% }-E^{**})+(\widetilde{I}-I^{**})\right]^{2}$

with Lyapunov derivative,

$\displaystyle\dot{\mathcal{F}}=\left[(\widetilde{S}-S^{**})+(\widetilde{E}-E^{% **})+(\widetilde{I}-I^{**})\right]\dot{\widetilde{N}}(t),$

where, $\widetilde{N}(t)=\widetilde{S}+\widetilde{E}+\widetilde{I}$ , it can be easily shown that $\dot{\widetilde{N}}(t)=\mu N_{0}-\mu\widetilde{N}(t)$ and $S^{**}+E^{**}+I^{**}=N_{0}$ , thus

$\displaystyle\dot{\mathcal{F}}=\left[(\widetilde{S}-S^{**})+(\widetilde{E}-E^{% **})+(\widetilde{I}-I^{**})\right]\left[\mu(S^{**}+E^{**}+I^{**})-\mu(% \widetilde{S}+\widetilde{E}+\widetilde{I})\right],$

so that

$\displaystyle\dot{\mathcal{F}}=-\mu\left[(\widetilde{S}-S^{**})+(\widetilde{E}% -E^{**})+(\widetilde{I}-I^{**})\right]^{2}.$ (56)

It follows from Eq. (56) that $\dot{\mathcal{F}}\leqslant$ 0 for $\mathcal{R}_{cons}>$ 1 with $\dot{\mathcal{F}}=$ 0 if and only if $\widetilde{S}=S^{**},\widetilde{E}=E^{**}$ and $\widetilde{I}=I^{**}$ . Thus, by LaSalleæŠ¯ Invariance Principle [11, 18] $(\widetilde{S},\widetilde{E},\widetilde{I})\rightarrow(S^{**},E^{**},I^{**})$ as $t\rightarrow\infty$ . Thus, the unique endemic equilibrium of the model Eq. (5) is GAS in $\mathcal{D}_{01}$ whenever $\mathcal{R}_{cons}>$ 1. ∎

6. Conclusion

An age structured of Herpes Simplex virus type 2 mathematical model is considered. A full stability analysis is performed. It showed that the model have a locally-asymptotically and a globally-asymptotically stable disease-free equilibrium point whenever the reproduction number $\mathcal{R}_{0}<$ 1. This result indicated that the spread of HSV-2 can be controlled whenever the initial sizes of the sub-populations of the model are in the basin of attraction of the DFE point. The existence of the endemic equilibrium point of the transformed model Eq. (2.1) is investigated. A unique endemic equilibrium point EEP of the model Eq. (2.1) is depicted whenever $\mathcal{R}_{0}>$ 1, this equilibrium is locally-asymptotically stable whenever $\mathcal{R}_{0}>$ 1. For the case where the effective contact rate and natural death rate are constants, it is shown that the model has a unique endemic equilibrium which is globally-asymptotically whenever the associated reproduction number exceeds unity.

Furthermore, it is shown that adding age-factor to the corresponding autonomous HSV2 model (considered in Podder and Gumel [23] for the case when quiescent infectious individuals can transmit infection at the same rate as the non-quiescent infectious individuals) does not alter the qualitative dynamics of the autonomous system (with respect to the elimination or persistence of the disease).

References

Balasubramaniam

Kuperstein

Stoopler

. Update on oral herpes virus infections. Dental Clinics of North America April 2014; 58(2): 265-280.

Blower

Porco

Darby

. Predicting and preventing the emergence of antiviral drug resistance in HSV-2. Nature Medicine 1998; 4: 673-678.

Chayavichitsilp

Buckwalter

Krakowski

Friedlander

. Herpes simplex. Pediatr Rev April 2009; 30(4): 119-129.

Elveback

, et al. Stochastic two-agent epidemic simulation models for a community of families. Amer J Epidemiol 1971; 267-280.

Fleming

McQuillan

Johnson

, et al. Herpes simplex virus type 2 in the United States, 1976 to 1994. N Engl J Med 1997; 337(16): 1105-11.

Gershengorn

Blower

. Impact of antivirals and emergence of drug resistance: HSV-2 epidemic control. AIDS Patient Care and STDs 2000; 14(3): 133142.

Greenhalgh

. Analytical results on the stability of age-structured epidemic models. IMA J Math Appl Med Biol 1987; 4: 109-144.

Greenhalgh

. Analytical thershold and stability results on age-structured epidemic models with vaccination. Theor Popul Biol 1987; 33: 266-290.

Gripenberg

. On a nonlinear integral equation modelling an epidemic in an age-structured. J Reine Angew Math 1983; 341: 54-67. population.

10.

Gurtin

MacCamy

. Non-linear age-dependent population dynamics. Arch Rational Mech Anal 1985; 54: 281-300.

11.

Hale

. Ordinary Differential Equations. New York: John Wiley and Sons; 1969.

12.

Hethcote

. The mathematics of infectious diseases. SIAM Rev 2000; 42: 599-653.

13.

Inaba

. Threshold and stability results for an age-structured epidemic model. J Math Biol 1990; 28: 411-434.

14.

Kamali

Nunn

Mulder

Van Dyck

Dobbins

Whitworth

. Seroprevalence and incidence of genital ulcer infection in a rural Ugandan population. Sex Transm Inf 1999; 75: 98-102.

15.

Katzmann

Dietz

. Evaluation of age-specific vaccination strategies. Theor Popul Biol 1984; 25: 125-137.

16.

Khan

Islam

Zaman

. Media coverage campaign in Hepatitis B transmission model. Appl Math Comp 2018; 331: 378-393.

17.

Agusto

Khan

. Optimal control strategies for dengue transmission in Pakistan. Math Biosc 2018; 305: 102-121.

18.

LaSalle

. The stability of dynamical systems. in: CBMS-NSF Regional Conf Ser in Appl Math, SIAM, Philadelphia. 1976.

19.

Levins

Awerbuch

Brinkman

Eckardt

Epstein

Makhoul

de Possas

Puccia

Spielman

Wilson

. The emergence of new diseases. American Scientist 1994; 82: 52-60.

20.

, et al. Thershold and stability results for an age-structured SEIR epidemic model. Comp Math Appl 2001; 42: 883-907.

21.

Liu

. Stability of an age-structured epidemiological model for hepatitis C. J Appl Math Comput 2008; 27: 159-173.

22.

Liu

Martcheva

. An age-structured two-strain epidemic model with super-infection. Math Biosci Engrg 2010; 7: 23-47.

23.

Podder

Gumel

. Qualitative dynamics of a vaccination model for HSV-2. IMA J Appl Math 2010; 75(1): 75-107.

24.

Royden

. Real Analysis. 3rd edition. Prentice Hall/Pearson Education; 1988.

25.

Safi

Gumel

Elbasha

. Qualitative analysis of an age-structured SEIR epidemic model with treatment. Appl Math Comput 2013; 219: 10627-10642.

26.

DarAssi

Safi

Al-Hdaibat

. A delayed SEIR epidemic model with pulse vaccination and treatment. Nonlinear Studies 2018; 3521-534.

27.

Safi

. Mathematical analysis of an age-structured quarantine/isolation model. Italian Journal of Pure and Applied Mathematics 2018; 39: 220-242.

28.

Schenzle

. An age-structured model of pre- and post-vaccination measles transmission. IMA J Math Appl Med Biol 1984; 1: 169-191.

29.

Smith

Waltman

. The Theory of the Chemostat. Cambridge University Press; 1995.

30.

Sutton

Banks

Castillo-Chavez

. Public vaccination policy using an age-structured model of pneumococcal infection dynamics. J Biol Dynamics 2009; 4(2): 176-195.

31.

Thieme

. Mathematics in Population Biology. Princeton: Princeton University Press; 2003.

32.

Obasi

Mosha

Quigley

, et al. Antibody to herpes simplex virus type 2 as a marker of sexual risk behaviour in rural Tanzania. J Infect Dis 1999; 179(1): 16-24.

33.

Wang

Zhang

Kuniya

. Global dynamics for a class of age-infection HIV models with nonlinear infection rate. Journal of Mathematical Analysis and Applications 2015; 432(1): 289-313.

34.

Wang

Zhang

Kuniya

. A note on dynamics of an age-of-infection cholera model. Mathematical Biosciences and Engineering 2016; 13(1): 227-247.

35.

Wang

Zhang

Kuniya

. The dynamics of an SVIR epidemiological model with infection age. IMA Journal of Applied Mathematics 2016; 81: 321-343.

36.

Wang

Lang

Zou

. Analysis of a structured HIV infection model with both virus-to-cell infection and cell-to-cell transmission. Nonlinear Analysis: Real World Applications 2017; 34: 75-96.

37.

Weiss

Buvé

Robinson

, et al. Study group on heterogeneity of HIV epidemics in African cities. The epidemiology of HSV-2 infection and its association with HIV infection in four urban African populations. AIDS 2001; 15(supp 4): S97-S108.

38.

Weiss

. Epidemiology of herpes simplex virus type 2 infection in the developing world. Herpes 2004; 11(Supplement 1).

39.

Webb

. Theory of Nonlinear Age-Dependent Population Dynamics. New York Basel: Dekker; 1985.