Transmission dynamics of two strain herpes simplex virus

1. Introduction

Herpes simplex virus (HSV) is a kind of infection that causes herpes. Herpes can appear in different part of the body, particularly on the mouth or genital area. HSV infections are endemic throughout the world[1, 2, 3, 4, 5]. For both point-prevalence and prospective studies, a large percentage of persons that are seropositive for HSV type 1 (HSV-1) or HSV type 2 (HSV-2) have no clinical manifestations of the disease [1, 2, 3, 4, 5,6]. HSV-1 and HSV-2 are lifelong infections [1, 2]. HSV-1 is mainly transmitted by oral to oral contact to cause infection in or around the mouth (oral herpes). HSV-2 is almost exclusively sexually transmitted, causing infection in the genital or anal area (genital herpes). Nonetheless, HSV-1 can also be transmitted to the genital area through oral-genital contact to cause genital herpes [2].

In \(2012\), it was estimated that about 3.7 billion people under the age of \(50\), or \(67\%\) of the population had HSV-1 infection while over \(267\) million women and \(150\) million men were living with HSV-2 infection [2]. Estimated prevalence of HSV-1 infection was highest in Africa, \(87\%\) and lowest in the Americas \(40-50\%\) [2]. Prevalence of HSV-2 infection was estimated to be highest also in Africa \((31.5\%)\), followed by the Americas \((14.4\%)\) [2, 7]. This is a clear sign that there should be a global call to fight against the global burden of HSV.

Many mathematical models have been developed to study the dynamics of HSV-2 see [8, 9, 10, 11, 12,13,14, 15]. In the manner that, Sally Blower and Li ma [13] formulated a mathematical model that predicts the effect of high prevalence of HSV-2 on HIV. Their results showed that HSV-2 epidemic has more than double impact on the peak of HIV incidence. Abu-Raddad et al [16], considered a homosexual male population and suggested that HSV-2 prevalence, if near endemic level may predict the spread of HIV. Foss et al. [17], developed a dynamical model to estimate the HIV infection due to HSV-2 in heterosexual population. Alvey et al. [18], developed a model that takes into account the transmission through either homosexual or heterosexual behaviour and investigate the impact of the coupled dynamics of HIV and HSV-2. Their results showed that homosexual transmission has great impact on the disease prevalence.

For all the models considered so far and to the best of our knowledge, none of them considered the transmission dynamics of the two strain of HSV. However, a detailed study of the transmission dynamics of disease with multiple strains has been one of the important problems in epidemiology. Hence, in order to fight against this global socio-economic burden, a deterministic model of two strains of HSV is developed and rigorously analysed, to get more insight into the long-term dynamics of the disease.

The paper is organized as follows. The model is formulated in Section 2 while it is analysed in Section 3. Also the numerical simulation is presented in Section 4 and the conclusion is drawn in Section 5.

2. Model formulation

The model is based on the transmission dynamics of HSV. The total population at time \(t\) is denoted by \(N(t)\). This is divided into susceptible individuals \(S(t)\); infectious individuals infected with strain \(1\) \(I_1(t)\), infectious individuals infected with strain \(2\)~ \(I_2(t)\); such that $$N(t)=S(t)+I_1(t)+I_2(t).$$ We give the detailed explanation of the transmission in the schematic diagram below:

It is assumed that the susceptible population is generated at a constant rate \(\wedge\). It is reduced due to infection by either HSV-1 (at a rate \(\beta_1\)) or by HSV-2 (at a rate \(\beta_2\)) and natural death at a rate \(\mu\), such that $$\frac{dS}{dt}=\wedge-\beta_1S I_1 -\beta_2 SI_2 -\mu S.$$ The rate of change of the infectious individuals infected with strain \(1\) is increased by infecting the susceptible population (at a rate \(\beta_1\)). It is decreased due to genital herpes induced by HSV-1 (at a rate \(\eta\)) and by natural death (at a rate \(\mu\)), so that $$\frac{dI_1}{dt}=\beta_1S I_1 -\left(\eta+\mu\right)I_1.$$ Finally, the rate of change of the infectious individuals infected with strain \(2\) is increased by infecting the susceptible population (at a rate \(\beta_2\)) and due to genital herpes induced by HSV-1 (at a rate \(\eta\)). It is decreased due natural death (at a rate \(\mu\)). Thus $$\frac{dI_2}{dt}=\beta_2 SI_1+ \eta I_1-\mu I_2.$$ In summary, the model for the transmission dynamics of both HSV-1 and HSV-2 is given by the non-linear differential equations. The model variables and parameters are described in 1 and the schematic diagram of the model is depicted in Figure 1.

Table 1. Description of model variables and parameters.

Variables	Description
\(S(t)\)	Susceptible Population
\(I_1(t)\)	Individuals infected with strain 1
\(I_2(t)\)	Individuals infected with strain 2
\(\wedge\)	Recruitment rate
\(\beta_1\)	Efficient contact rate of strain 1
\(\beta_2\)	Efficient contact rate of strain 2
\(\eta\)	Rate of genital herpes induced by strain 1
\(\mu\)	Natural mortality rate

2.1. Basic properties of the model

For the model to be epidemiologically meaningful, we need to prove that all the states variables are non-negative for all time \(t\). In other words the solution of Equation (1) with positive initial condition will remain positive for all \(t\geq 0\).

Theorem 1. If the initial condition \(S(0)>0\), \(I_1(0)>0\) and \(I_2(0)>0\), then the solutions \(S\), \(I_1\), and \(I_2\) of model Equation (1) are positive for all \(t\).

Proof. Given that \(S(0)>0\), \(I_1(0)>0\) and \(I_2(0)>0\), we want to show that \(S(t)>0\), \(I_1(t)>0\) and \(I_2(t)>0\). Let \(\tau=sup\left\{t>0: S(t)>0, I_1(t)>0, I_2(t)>0\right \}\). Thus \(\tau >0\). Lets consider the first equation in Equation (1), we have that

Integrating the above equation, we get Considering the second equation in Equation (1), we have that Integrating the above equation, we get Finally, considering the third equation in Equation (1), we have that Integrating the above equation, yields Thus \(S(t)\), \(I_1(t)\) and \(I_2(t)\) are all positive for any non-negative initial conditions.

2.2. Positively invariant region

Lemma 2. The region $$\Gamma=\left\{\left( S,I_1,I_2\right) \in \mathbb{R}^{3}_{+}: S+I_1+I_2\leq \frac{\wedge}{\mu}\right \} $$ is positively invariant and attracting in Equation (1).

Proof. Adding the Equation (1), gives \begin{align*} \frac{dN}{dt}&=\wedge-\mu\left(S+I_1+I_2\right),\\ \frac{dN}{dt}&\leq\wedge-\mu N. \end{align*} Thus the population is bounded above by \(\frac{\wedge}{\mu}\) so that \(\frac{dN}{dt}\frac{\wedge}{\mu}\). Hence by comparison Theorem [19], \(N(t)\leq N(0)e^{-\mu t}+\frac{\wedge}{\mu}\left[1-e^{-\mu t}\right]\), in particular \(N(t)\leq \frac{\wedge}{\mu}\) if \(N(0)\leq \frac{\wedge}{\mu}\). Thus, \(\Gamma\) is positively invariant. Also if \(N(0)>\frac{\wedge}{\mu}\), then either the solution enters \(\Gamma\) in finite time or \(N(t)\) approaches \(\frac{\wedge}{\mu}\) asymptotically. Therefore, the region \(\Gamma \) attracts all solution in \(\mathbb{R}^{3}_{+}\). Hence, it suffices to consider the dynamics of the model Equation (1) in \(\Gamma\), where the model is epidemiologically and mathematically well posed [20].

3. Model analysis

3.1. Disease-free equilibrium point

At equilibrium point, we equate each of the right hand side of Equation (1) to zero and solve. Where we represent diseases free equilibrium (DFE) as \(E_0\) to be

3.2. Basic reproduction number (\(\mathcal{R}_0\))

The basic reproduction number \(\mathcal{R}_0\) measures the average number of secondary new infections caused by one primary infection in a completely susceptible population. \(\mathcal{R}_0\) gives the threshold whether a disease will go extinction or persist. In this research, we find our \(\mathcal{R}_0\) using the next generation matrix method on Equation (1). Using the notation in [21], the matrices \(\mathscr{F}\) is the rate of appearance of new infections in compartment \(i\), and \(\mathscr{V}\) is the rate of other transitions between compartment \(i\) and other infected compartments of Equation (1) are given respectively by \begin{equation*} \mathscr{F}= \begin{bmatrix} \beta_1SI_1\\ \beta_2SI_2+\eta I_1 \end{bmatrix}, \end{equation*} \begin{equation*} \mathscr{V}= \begin{bmatrix} \left(\eta +\mu\right)\mu\\ \mu I_2 \end{bmatrix}. \end{equation*} Computing the matrices \(F\) and \(V,\) for the new infection terms and of the transition terms, respectively, we have \begin{equation*} F= \begin{bmatrix} \frac{\beta_1 \wedge}{\mu}&0\\ \eta& \frac{\beta_2 \wedge}{\mu} \end{bmatrix}, \end{equation*} \begin{equation*} V= \begin{bmatrix} \eta +\mu&0\\ 0& \mu \end{bmatrix}. \end{equation*} Thus, the basic reproduction number of Equation (1), denoted by \(\mathcal{R}_0\), is given by (where \(\rho\) is the spectral radius) \begin{equation*} \mathcal{R}_0=\rho\left(FV^{-1}\right)=max\left\{\mathcal{R}_{1},\mathcal{R}_{2}\right\}, \end{equation*} where \(\mathcal{R}_{1}\) and \(\mathcal{R}_{2}\) are the associated reproduction numbers of strain \(1\) and strain \(2\), respectively given by

3.3. Local stability of disease-free equilibrium

Theorem 3. The DFE of Equation (1), given by \(E_0\), is locally asymptotically stable (LAS) if \(\mathcal{R}_01\).

Proof. The Jacobian matrix of Equation (1) denoted by \(J\) is given as

evaluating \(J\) at \(E_0\), we get \begin{equation*} J_{\left(\frac{\wedge}{\mu},0,0\right)}= \begin{bmatrix} -\mu&-\frac{\beta_1 \wedge}{\mu}&-\frac{\beta_2 \wedge}{\mu}\\ 0& \frac{\beta_1 \wedge}{\mu}-\left(\eta +\mu\right)&0\\ 0& \eta& \frac{\beta_2 \wedge}{\mu}-\mu \end{bmatrix}. \end{equation*} The eigenvalues associated with the above matrix are \begin{align*} \lambda_1&= -\mu,\\ \lambda_2&=\frac{1}{\mu}\left(\beta_1 \wedge-\eta \mu-\mu^2\right),\\ \lambda_3&= \frac{1}{\mu}\left(\beta_2 \wedge -\mu^2\right). \end{align*} But \(\lambda_2\) and \(\lambda_3\) can be re-written as \(\left(\wedge+\mu\right)\left(\mathcal{R}_{1}-1\right)\) and \(\frac{1}{\mu}\left(\mathcal{R}_{2}-1\right)\) respectively. It is clear that since \(\lambda_1< 0\), then \(E_0\) is LAS if \(\lambda_2< 0\) and \(\lambda_3< 0\) i.e if \(\mathcal{R}_{1}< 1\) and \(\mathcal{R}_{2}1\),\(\mathcal{R}_{2}>1\) or both \(\mathcal{R}_{1}\), \(\mathcal{R}_{2}\) \(>0\).

3.4. Global stability of disease-free equilibrium

Theorem 4. The DFE given by \(E_0\) is globally asymptotically stable (GAS) in \(\Gamma\) whenever \(\mathcal{R}_{0}\leq 1\).

Proof. Consider the Lyapunov function for Equation (1)

Clearly the function \(\mathcal{F}\) is continuous and positive definite, with lyapunov derivative given below $$\dot{\mathcal{F}}=K_1\dot{I_1}+K_2\dot{I_2}.$$ Substituting \(\dot{I_1}\) and \(\dot{I_2}\) in the above equation, we get where $$K_1=\left(\mu \eta +\beta_1 \beta_2\right)^4 ~~and~~K_2=\left(\frac{\beta_1 \wedge}{\mu+\eta}-\mu\right)^2.$$ Thus from Equation (13), \(\dot{\mathcal{F}}\leq 0\) if \(\mathcal{R}_{0}={max}\left\{\mathcal{R}_{1},\mathcal{R}_{2}\right\}\leq 1\) and \(\dot{\mathcal{F}}=0\) if and only if \(I_1=I_2=0.\) Substituting \(I_1=I_2=0\) in Equation (1) we get \(S(t)\to \frac{\wedge}{\mu} ~~{as}~~ t\to \infty\). Also the largest compact invariant set in \(\left\{\left( S,I_1,I_2\right)\in \Gamma: \frac{d\mathcal{F}}{dt}=0\right\}\) is the singleton \(\left\{E_0\right\}.\) Hence it follows from the Laselle invariance principle [22], every solutions in Equation (1) with initial condition in \(\mathbb{R}^{3}_{+}\), converges to DFE \(E_0\) as \(t\to \infty\) whenever \(\mathcal{R}_{0}=\text{max}\left\{\mathcal{R}_{1},\mathcal{R}_{2}\right\}\leq 1\). Thus the DFE, \(E_0\), is GAS in \(\Gamma\) for \(\mathcal{R}_{0}=\text{max}\left\{\mathcal{R}_{1},\mathcal{R}_{2}\right\}\leq 1\).

3.4.1. Existence of endemic equilibrium

Finding the condition(s) for the existence of an equilibrium \(E_* =\left(S^{*},I^{*}_1, I^{*}_2\right)\) of the model Equation (1) such that the disease is endemic in the population, the Equation (1) are solved at steady state which yields

3.5. Local stability of endemic equilibrium

Theorem 5. The EE, \(E_1\) of Equation (1) is LAS provided that

where \(q_1\), \(q_2\), \(q_3\), \(q_4\), \(q_5\), and \(q_5\) are given in the proof.

Proof. The characteristics polynomial of Equation (10) at \(E_1\) can be written as

where with Under the conditions in Equation (14), (15), (16) and (17), \(a_1>0\), \(a_2>0\) and \(a_3>0\). But \(a_1 a_2=q_1q_3+q_2q_4-q_1q_4-q_2q_3\) and \(a_1 a_2>a_3\) if \(q_1q_3+q_2q_4-q_1q_4-q_2q_3>q_5-q_6\) that is if Thus by the Routh- Hurwitz criterion [23] all the eigenvalues of Equation (10) at \(E_1\) have negative real part. Hence the endemic equilibrium of Equation (1) subject to the conditions of Equation (14), (15), (16) ,(17) and (20) is LAS.

3.6. Global stability of endemic equilibrium: special case

Theorem 6. If \(\mathcal{R}_{0}=\text{max}\left\{\mathcal{R}_{1},\mathcal{R}_{2}\right\}\geq1\) and \(\eta=0\), then the endemic equilibrium is globally asymptotically stable (GAS) in \(\Gamma\).

Proof. Consider the non-linear Lyapunov function (such non-linear functions has been used in many mathematical epidemiology, see [24, 25, 26, 27] for

with Lyapunov derivative $$\dot{\mathcal{F}}=\dot{S}-\frac{S^{*}}{S}\dot{S}+\dot{I_1}-\frac{I^{*}_1}{I_1}\dot{I_1}+\dot{I_2}-\frac{I^{*}_2}{I_2}\dot{I_2}.$$ Substituting \(\dot{S}\), \(\dot{I_1}\) and \(\dot{I_2}\) in the above equation, we get It can be shown from Equation (1) with \(\eta=0\), at endemic state that \begin{align}\label{w1} \wedge&=\beta_1I^{*}_1S^*+\beta_2I^{*}_2S^*+\mu S^*, %\mu I^{*}_1&=\beta_1I^{*}_1S^*\\ %\mu I^{*}_2&=\beta_2I^{*}_2S^* \end{align} substituting Equation (23) into (22), we have It can be shown from Equation (1) with \(\eta=0\), at endemic state that substituting Equation (25) into (24), we have from Equation (25), we have (26) to be Since the arithmetic mean is greater than or equal to the geometric mean, then \begin{equation*} \begin{split} 2-\frac{S^*}{S}-\frac{S}{S^*}&\leq 0. \end{split} \end{equation*} Clearly, the first term of \(\dot{\mathcal{F}}\) is always negative and the second and third terms are also negative, therefore \(\dot{\mathcal{F}}\leq 0,\) with \(\dot{\mathcal{F}}=0\) if and only if \(S^*=S\). Also the largest compact invariant set in \(\left\{\left( S,I_1,I_2\right)\in \Gamma: \frac{d\mathcal{F}}{dt}=0\right\}\) is the singleton \(\left\{E_1\right\}.\) Thus by the Laselle Invariance Principle [22], every solutions in Equation (1) with initial conditions in \(\mathbb{R}^{3}_{+}\), converges to \(E_1\) as \(t\to \infty\) whenever \(\mathcal{R}_{0}=max\left\{\mathcal{R}_{1},\mathcal{R}_{2}\right\}\geq 1\). Hence, \(E_1\), is GAS in \(\Gamma\) for \(\mathcal{R}_{0}=max\left\{\mathcal{R}_{1},\mathcal{R}_{2}\right\}> 1\).

4. Numerical simulations and discussions

In this section, the model Equation (1), using the parameter values given in Table 2 (unless otherwise stated), to assess the impact of super-infection on the dynamics of two strain HSV. The objective of this section is to illustrate some of the theoretical results in this paper. Since the model presented in this paper is completely new (no similar two strain model for HSV model has yet been published in the literature to our knowledge), appropriate data for estimating the associated parameters are not available at the present time. Thus, the parameter values chosen for the numerical simulations below may not all be realistic biologically, although such uncertainties in parameter values are partially addressed below by considering different rate of genital herpes induced by strain 1(super-infection) in the simulations, it is important to emphasize that the simulation results obtained should be interpreted bearing these uncertainties in mind.

Table 2. Parameter values for model Equation (1).

Parameters	Baseline value	Reference
\(\wedge\)	10000 per day	Assumed
\(\beta_1\)	\(7\times 10^{-9}\) per day	Assumed
\(\beta_2\)	\(2\times 10^{-9}\) per day	Assumed
\(\eta\)	\(0.0022\)  per day	Assumed
\(\mu\)	\(0.0167\)  per day	Assumed

With the given parameters in Table 2, \(\mathcal{R}_{1}=0.204\) and \(\mathcal{R}_{2}=0.072\) such that \(\mathcal{R}_{0}=0.204< 1\). It is observed that the susceptible individuals increases as a result \(\wedge\) which is the recruitment rate, however due to natural mortality rate \(\mu\) the increase in susceptible individuals is regulated, moreover individuals in the susceptible compartment reach a saturation point \(600000\) (\(\frac{\wedge}{\mu}\)) from Figure (2a). This shows that as time increases the individuals who enter the susceptible compartment would equal to those who leave the compartment as a result of natural mortality rate \(\mu\). It was also observed that whenever there is no disease Figure (2b) and (2c) goes to zero. Therefore by Theorem 4 the DFE is GAS. Figure (2a), (2b) and (2c) shows this this simulations, confirming the GAS property of the DFE. In addition, the effect of super-infection governed by \(\eta\) is monitored, which shows that the total number of infected individuals with strain \(2\) increases with increase in \(\eta\) and decreases with decrease in \(\eta\). This is depicted in Figure (2b) and (2c). Furthermore, in Figure (2d) and (3), additional simulation shows that when \(\mathcal{R}_{2}< 1< \mathcal{R}_{1}\), strain \(1\) does not drive strain \(2\) to extinction but when \(\mathcal{R}_{1}< 1< \mathcal{R}_{2}\), strain \(2\) drives out strain \(1\) to extinction (competitive exclusion).

Figure (4a) depicts the scenario where the susceptible individuals increases then decreases but stabilizes in the system as a result of \(\mathcal{R}_{0}>1\) and Figure (4b) and (4c) shows how individuals infected with srain \(1\) and strain respectively evolve due to \(\mathcal{R}_{0}>1\). Also, Figure (4b) and (4c) shows that when \(\mathcal{R}_{1}>1\) and \(\mathcal{R}_{2}>1\), strain \(1\) and strain 2 co-exist.

Conclusions

A new deterministic model for the transmission dynamics of two-strain HSV is designed and analysed. The main findings in this paper are:

1. Model Equation (1) has a GAS DFE whenever \(\mathcal{R}_0< 1.\)
2. Model Equation (1) has a GAS EE for special case when \(\eta=0\) whenever \(\mathcal{R}_0>1.\)
3. Numerical simulation of model Equation (1) shows that strain 2 drives strain 1 to extinction when \(\mathcal{R}_1< 1< \mathcal{R}_2\), that is the model undergoes competitive exclusion. But when \(\mathcal{R}_2< 1< \mathcal{R}_1\), strain 1 those not drive strain 2 to extinction.
4. Numerical simulation of Equation (1) shows that super-infection has a great impact on strain 2 by increasing it population for any small increment on the super-infection parameter \(\eta.\)
5. Numerical simulation of Equation (1) shows that the two strains coexist when \(\mathcal{R}_1>\mathcal{R}_2>1\).

However, this study shows that HSV-1 has great impact on increasing the spread of HSV-2. This is as the result of the genital herpes induced by HSV-1 (super-infection). This study suggest that governments should provide screening centres to enable every individual to know his or her status of HSV. Also the governments and public health agencies should organize massive campaign world wide to sensitize people on HSV and to encourage the general public to go for HSV testing as this would enable researchers get ample data to carry out research effectively and provide possible solutions to eradicating this burden (HSV).

Acknowledgments

One of the authors Shamsuddeen Ibrahim acknowledges, with thanks the support of the MasterCard Foundation, the Canadian Government and African Institute for Mathematical Sciences Ghana. The authors are grateful to the anonymous Reviewers and the Handling Editor for their constructive comments.

Author Contributions

All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.

Competing Interests

The author(s) do not have any competing interests in the manuscript.