Modeling the Transmission Dynamics of the Monkeypox Virus Infection with Treatment and Vaccination Interventions

Show more

1. Introduction

Monkeypox is an infectious disease caused by the monkeypox virus [1] . Monkeypox virus is a zoonotic viral disease that occurs primarily in remote villages of Central and West Africa in proximity to Tropical Rainforest where there are more frequent contacts with infected animals [2] . The monkeypox virus which is closely related to variola virus was first identified by Magnus et al. (1959) as the causal agent in two outbreaks of pox infections in cynomolgus monkeys that were then received by the Statens Serum-institut, Copenhagen, Denmark, from Singapore [3] [4] . Human monkeypox is clinically related to smallpox as the two infections are difficult to distinguish from [2] . Monkeypox is usually transmitted to humans from rodents, pets and primates through contact with animal blood or bites [2] [5] . The infection can also be found in Gambian pouched rats (Cricetomys gambianus), dormice (Graphiurus sp.) and African squirrels (Heliosciurus and Funisciurus) [2] . Transmission of monkeypox virus occurs when a person comes into contact with the virus from an infected animal, human or material contaminated with the virus [5] . The virus enters the body through broken skin (even if not visible), respiratory tract or the mucous membranes (eyes, nose or mouth) [5] . Animal to human transmission may occur by bite or scratch, bush meat preparations direct contact with body fluids or lesion material such as through contaminated bedding. Human to human transmission is thought to occur primarily through large respiratory droplets which generally cannot travel more than few feet, and therefore, prolonged face-to-face contact is required [5] . Other methods of transmission in this category include direct contact with body fluids or lesion material materials [5] .

In humans, the symptoms of the monkeypox virus infection are similar to but milder than the symptoms of smallpox. The infection usually begins with fever, headache, muscle ache and exhaustion [5] [6] . The main difference between symptoms of smallpox and monkeypox is that monkeypox causes lymph nodes to swell (lymphadenopathy) while smallpox does not [5] . The incubation period for monkeypox is usually 7 - 14 days but can range from 5 - 21 days [5] .

Evidence of viral infection in humans with monkeypox virus was first identified in the Democratic Republic of Congo (DRC), formerly known as Zaire, in the town of Basankusu, Equateur Province in the year 1970 [3] . A second outbreak of the infection occurred in the DRC/Zaire in 1996-1997 [5] . In 2003, a small outbreak of human monkeypox virus in United States occurred among owners of pet prairie dogs [5] . The outbreak originated from Villa Park, Illinois, outside of Chicago, when an exotic animal dealer kept young prairie dogs in close proximity to an infected Gambian pouched rat imported from Ghana, where a total of 71 people were reportedly infected with no deaths [5] . A more recently, an outbreak of the monkeypox virus infection has occurred in the Bayelsa State of the Federal Republic of Nigeria [6] . The Nigeria Centre for Disease Control (NCDC) has confirmed 4 cases out of the 43 suspected samples of the monkeypox virus infection sent the World Health Organization (WHO) Regional Laboratory in Dakar, Senegal, with 32 close contacts of the cases placed under clinical watch [6] [7] .

Between 1970 and 1986, 404 cases of monkeypox virus were reported in 7 West and Central Africa countries (DRC, Ivory Coast, Sierra Leone, Cameroon, Central African Republic, Liberia and Nigeria) among mainly children [8] . There are several promising antiviral drugs under development which may offer therapeutic benefit for monkeypox patients of which Cidofovir has demonstrated protection in challenge studies performed under animal models [8] . There is no known monkeypox vaccine in circulation as at now, but the smallpox vaccine (vaccinia) which has demonstrated protection against monkeypox with about 85% vaccine efficacy could be used in monkeypox-endemic areas [8] . The known complications of the smallpox vaccine, that is, the increase of the HIV/AIDS prevalence in monkeypox-endemic environments are to be watched when using the smallpox vaccine for immunization [8] .

Currently, there is not much on the modeling aspect of the monkeypox virus infection [9] . But, the work in [9] provides framework for studying the transmission dynamics of the pox-like viruses with monkeypox as case study where they divided the host into primates transmitting the virus to humans through contact with infected rodents, which is a probability function, and human-to-human with increased transmission rate which is also a probability function of contacts with infected rodents or humans. Another mathematical modeling work for the transmission dynamics of smallpox virus with control interventions can be traced to [10] , where they studied the dynamics of the virus in human host only, sustaining the virus from contact with an infected rodent or human. They studied the dynamics on different modeling schemes of SIR and SEIR approaches. In [11] , control measures on respiratory pathogens may include any or all of the policies; quarantine, infection control precautions, case identification and isolation and immunization interventions.

Therefore, this paper is set out to review the existing work of [9] by incorporating control interventions of treatment and vaccination, and latency/exposure period on the trends of successive chains of progression in both primates and human hosts since the monkeypox virus infection has incubation rates in the humans [5] and the primates (on assumption).

2. Model Formulations

2.1. Description of the Model

The model in this paper divides the host population into two; the non-human primates and/or some wild rodents, and the humans host population. The non-human primate population was further divided into Susceptible $\left({S}_{n}\right)$ , Exposed/Latent $\left({E}_{n}\right)$ , Infected $\left({I}_{n}\right)$ and Recovered $\left({R}_{n}\right)$ subpopulations. The non-human primates and/or some wild animals are recruited into the Susceptible $\left({S}_{n}\right)$ class at a constant birth rate ${\Lambda}_{n}$ and become exposed to the monkeypox virus after getting into contact with an infected non-human primate at a rate ${\lambda}_{n}$ with

${\lambda}_{n}={\beta}_{n1}\frac{{I}_{n}}{{N}_{n}}$ (1)

where ${\beta}_{n1}$ is the product of effective contact rate and probability of the non-human primate getting infected per contact. After incubation of the virus is achieved, the exposed primate proceeds to the infected class $\left({I}_{n}\right)$ at a rate ${\nu}_{n}$ . The infected animals in $\left({I}_{n}\right)$ are capable of either; infecting other animals when they come into contact, die due to the disease at a rate ${d}_{n}$ or recover naturally with permanent immunity at a rate ${\rho}_{n}$ and move into $\left({R}_{n}\right)$ . All non-human primates in the model experience natural mortality rate ${\mu}_{n}$ .

The total human host population was also divided into Susceptible $\left({S}_{h}\right)$ , Vaccinated $\left({V}_{h}\right)$ , Exposed/Latent $\left({E}_{h}\right)$ , Infected $\left({I}_{h}\right)$ and Recovered $\left({R}_{h}\right)$ human subpopulations. The Susceptible humans are recruited into $\left({S}_{h}\right)$ through birth and migration at a constant rate ${\Lambda}_{h}$ . A susceptible individual is either vaccinated against the monkeypox virus at a rate ${\alpha}_{h}$ and move to $\left({V}_{h}\right)$ with permanent immunity or become exposed to the monkeypox virus after getting into contact with an infected human or non-human primates at a rate ${\lambda}_{h}$ with

${\lambda}_{h}={\beta}_{n2}\frac{{I}_{n}}{{N}_{n}}+{\beta}_{h}\frac{{I}_{h}}{{N}_{h}}$ (2)

where ${\beta}_{n2}$ is the product of the effective contact rate and probability of the human being infected per contact with an infectious non-human primate animal, and ${\beta}_{h}$ is the product of the effective contact rate and the probability of the human being infected with monkeypox virus after getting into contact with an infectious human per contact. After the incubation period, the Exposed human in $\left({E}_{h}\right)$ proceeds to the infected class $\left({I}_{h}\right)$ at a rate ${\nu}_{h}$ . Individuals in $\left({I}_{h}\right)$ either die due to the virus at a constant rate ${d}_{h}$ or recover with permanent immunity after receiving treatment at a rate ${\rho}_{h}$ into $\left({R}_{h}\right)$ . All individuals in the human subpopulations suffer a natural mortality at a constant rate ${\mu}_{h}$ . All parameters in the model are strictly nonnegative and will assume values presented in Table 1 during simulations and sensitivity analysis. The schematic diagram of the model in this description is presented in Figure 1 below.

2.2. Model Equations

From the description of the model and the schematic diagram presented in Figure 1 above, we derived the following model equations

${{S}^{\prime}}_{n}={\Lambda}_{n}-\left({\mu}_{n}+{\lambda}_{n}\right){S}_{n}$ (3)

${{E}^{\prime}}_{n}={\lambda}_{n}{S}_{n}-\left({\mu}_{n}+{\nu}_{n}\right){E}_{n}$ (4)

${{I}^{\prime}}_{n}={\nu}_{n}{E}_{n}-\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right){I}_{n}$ (5)

Table 1. Model parameter values.

Figure 1. Schematic model diagram.

${{R}^{\prime}}_{n}={\rho}_{n}{I}_{n}-{\mu}_{n}{R}_{n}$ (6)

${{S}^{\prime}}_{h}={\Lambda}_{h}-\left({\mu}_{h}+{\lambda}_{h}+{\alpha}_{h}\right){S}_{h}$ (7)

${{V}^{\prime}}_{h}={\alpha}_{h}{S}_{h}-{\mu}_{h}{V}_{h}$ (8)

${{E}^{\prime}}_{h}={\lambda}_{h}{S}_{h}-\left({\mu}_{h}+{\nu}_{h}\right){E}_{h}$ (9)

${{I}^{\prime}}_{h}={\nu}_{h}{E}_{h}-\left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right){I}_{h}$ (10)

${{R}^{\prime}}_{n}={\rho}_{h}{I}_{h}-{\mu}_{h}{R}_{h}$ (11)

${N}_{n}\left(t\right)={S}_{n}+{E}_{n}+{I}_{n}+{R}_{n}$ (12)

${N}_{h}\left(t\right)={S}_{h}+{V}_{h}+{E}_{h}+{I}_{h}+{R}_{h}$ (13)

Subject to the following nonnegative initial conditions:

${S}_{n}\left(0\right)\ge 0,{E}_{n}\left(0\right)\ge 0,\text{\hspace{0.17em}}{I}_{n}\left(0\right)\ge 0,\text{\hspace{0.17em}}{R}_{n}\left(0\right)\ge 0$ (14)

${S}_{h}\left(0\right)\ge 0,{V}_{h}\left(0\right)\ge 0,{E}_{h}\left(0\right)\ge 0,\text{\hspace{0.17em}}{I}_{h}\left(0\right)\ge 0,\text{\hspace{0.17em}}{R}_{h}\left(0\right)\ge 0$ (15)

$\begin{array}{l}{S}_{n}\left(0\right)+{E}_{n}\left(0\right)+{I}_{n}\left(0\right)+{R}_{n}\left(0\right)\le {N}_{n}\left(0\right)\text{\hspace{0.17em}}\text{and}\\ {S}_{h}\left(0\right)+{V}_{h}\left(0\right)+{E}_{h}\left(0\right)+\text{\hspace{0.17em}}{I}_{h}\left(0\right)+\text{\hspace{0.17em}}{R}_{h}\left(0\right)\le {N}_{h}\left(0\right)\end{array}$ (16)

3. Model Analysis

The model analysis begins by showing that all feasible solutions of the model are uniformly bounded in a proper subset of $\Omega $ . Thus the feasible region

$\Omega =\{\begin{array}{l}\left({S}_{n},{E}_{n},{I}_{n},{R}_{n}\right)\in {\mathbb{R}}_{+}^{4}\text{\hspace{0.17em}}:\text{\hspace{0.17em}}{N}_{n}\le \frac{{\Lambda}_{n}}{{\mu}_{n}}\\ \left({S}_{h},{V}_{h},{E}_{h},{I}_{h},{R}_{h}\right)\in {\mathbb{R}}_{+}^{5}\text{\hspace{0.17em}}:\text{\hspace{0.17em}}{N}_{h}\le \frac{{\Lambda}_{h}}{{\mu}_{h}}\end{array}$ (17)

is considered. Therefore, after differentiation of (12) and (13), and proper substitutions, we have:

$\frac{\text{d}{N}_{n}\left(t\right)}{\text{d}t}={\Lambda}_{n}-{\mu}_{n}{N}_{n}-{d}_{n}{I}_{n}\le {\Lambda}_{n}-{\mu}_{n}{N}_{n}$ (18)

and;

$\frac{\text{d}{N}_{h}\left(t\right)}{\text{d}t}={\Lambda}_{h}-{\mu}_{h}{N}_{h}-{d}_{h}{I}_{h}\le {\Lambda}_{h}-{\mu}_{h}{N}_{h}$ (19)

Applying [12] on the differential inequalities in (18) and (19), we obtained:

$\{\begin{array}{l}{N}_{n}\left(t\right)\le {N}_{n}\left(0\right){\text{e}}^{-{\mu}_{n}t}+\frac{{\Lambda}_{n}}{{\mu}_{n}}\left(1-{\text{e}}^{-{\mu}_{n}t}\right)\\ {N}_{h}\left(t\right)\le {N}_{h}\left(0\right){\text{e}}^{-{\mu}_{h}t}+\frac{{\Lambda}_{h}}{{\mu}_{h}}\left(1-{\text{e}}^{-{\mu}_{h}t}\right)\end{array}$ (20)

where ${N}_{n}\left(0\right)$ and ${N}_{h}\left(0\right)$ are the initial populations of the non-human

primates and the humans respectively. Therefore, $0\le {N}_{n}\le \frac{{\Lambda}_{n}}{{\mu}_{n}}$ and $0\le {N}_{h}\le \frac{{\Lambda}_{h}}{{\mu}_{h}}$ as $t\to \infty $ . This implies that, $\frac{{\Lambda}_{n}}{{\mu}_{n}}$ and $\frac{{\Lambda}_{h}}{{\mu}_{h}}$ are upper bounds for ${N}_{n}\left(t\right)$ and ${N}_{h}\left(t\right)$ respectively, as long as ${N}_{n}\left(0\right)\le \frac{{\Lambda}_{n}}{{\mu}_{n}}$ and ${N}_{h}\left(0\right)\le \frac{{\Lambda}_{h}}{{\mu}_{h}}$ . Hence, the

feasible solution of the model equations in (3)-(13) enters the region $\Omega $ which is a positively invariant set. Thus, the system is mathematically and epidemiologically well-posed. Therefore, for an initial starting point $x\in \Omega $ , the trajectory lies in $\Omega $ , and so it is sufficient to restrict our analysis on $\Omega $ . Clearly, under the dynamics described by the model equations, the closed set $\Omega $ is hence a positively invariant set.

3.1. Model Equilibrium Points

Using standard approaches, the model disease-free ${\epsilon}_{0}$ and endemic ${\epsilon}_{\ast}$ (which existed only when ${R}_{0}>1$ ) equilibrium points are established as follows:

$\begin{array}{c}{\epsilon}_{0}=\left({S}_{h},{V}_{h},{E}_{h},{I}_{h},{R}_{h},{S}_{n},{E}_{n},{I}_{n},{R}_{n}\right)\\ =\left(\frac{{\Lambda}_{h}}{{\alpha}_{h}+{\mu}_{h}},\frac{{\Lambda}_{h}}{{\mu}_{h}}\frac{{\alpha}_{h}}{{\alpha}_{h}+{\mu}_{h}},0,0,0,\frac{{\Lambda}_{n}}{{\mu}_{n}},0,0,0\right)\end{array}$ (21)

${\epsilon}_{\ast}=\left({S}_{h}^{\ast},{V}_{h}^{\ast},{E}_{h}^{\ast},{I}_{h}^{\ast},{R}_{h}^{\ast},{S}_{n}^{\ast},{E}_{n}^{\ast},{I}_{n}^{\ast},{R}_{n}^{\ast}\right)$ (22)

where:

${S}_{h}^{\ast}=\frac{{\Lambda}_{h}}{\left({\mu}_{h}+{\alpha}_{h}+{\lambda}_{h}^{\ast}\right)}$ , ${V}_{h}^{\ast}=\frac{{\Lambda}_{h}}{{\mu}_{h}}\frac{{\alpha}_{h}}{\left({\mu}_{h}+{\alpha}_{h}+{\lambda}_{h}^{\ast}\right)}$ , ${E}_{h}^{\ast}=\frac{{\lambda}_{h}^{\ast}{\Lambda}_{h}}{\left({\mu}_{h}+{\nu}_{h}\right)\left({\mu}_{h}+{\alpha}_{h}+{\lambda}_{h}^{\ast}\right)}$ ,

${I}_{h}^{\ast}=\frac{{\nu}_{h}{\lambda}_{h}^{\ast}{\Lambda}_{h}}{\left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right)\left({\mu}_{h}+{\nu}_{h}\right)\left({\mu}_{h}+{\alpha}_{h}+{\lambda}_{h}^{\ast}\right)}$ ,

${R}_{h}^{\ast}=\frac{{\Lambda}_{h}}{{\mu}_{h}}\frac{{\rho}_{h}{\nu}_{h}{\lambda}_{h}^{\ast}}{\left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right)\left({\mu}_{h}+{\nu}_{h}\right)\left({\mu}_{h}+{\alpha}_{h}+{\lambda}_{h}^{\ast}\right)}$ ,

${S}_{n}^{\ast}=\frac{{\Lambda}_{n}}{\left({\mu}_{n}+{\lambda}_{n}^{\ast}\right)}$ , ${E}_{n}^{\ast}=\frac{{\lambda}_{n}^{\ast}{\Lambda}_{n}}{\left({\mu}_{n}+{\nu}_{n}\right)\left({\mu}_{n}+{\lambda}_{n}^{\ast}\right)}$ , ${I}_{n}^{\ast}=\frac{{\nu}_{n}{\lambda}_{n}^{\ast}{\Lambda}_{n}}{\left({\mu}_{n}+{\nu}_{n}\right)\left({\mu}_{n}+{\lambda}_{n}^{\ast}\right)\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)}$ ,

${R}_{n}^{\ast}=\frac{{\Lambda}_{n}}{{\mu}_{n}}\frac{{\rho}_{n}{\nu}_{n}{\lambda}_{n}^{\ast}}{\left({\mu}_{n}+{\nu}_{n}\right)\left({\mu}_{n}+{\lambda}_{n}^{\ast}\right)\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)}$ with

${\lambda}_{n}^{\ast}={\beta}_{n1}\frac{{I}_{n}^{\ast}}{{N}_{n}^{\ast}}$ , ${\lambda}_{h}^{\ast}={\beta}_{n2}\frac{{I}_{{}_{n}}^{\ast}}{{N}_{n}^{\ast}}+{\beta}_{h}\frac{{I}_{{}_{h}}^{\ast}}{{N}_{h}^{\ast}}$ ,

${N}_{n}^{\ast}=\frac{{\Lambda}_{n}-{d}_{n}{I}_{n}^{\ast}}{{\mu}_{n}}$ and ${N}_{h}^{\ast}=\frac{{\Lambda}_{h}-{d}_{h}{I}_{h}^{\ast}}{{\mu}_{h}}$ .

3.2. Local Stability Analysis of the Model

3.2.1. Local Stability of the Disease-Free Equilibrium (DFE) Point

The basic reproduction number of the model was computed using the next-generation matrix as defined in [13] and [14] . It is defined to be largest eigenvalue or spectral radius of the characteristic equation $\left|F{V}^{-1}-\psi I\right|=0$ . Using the notations in [13] for the model system (3)-(11), the associated matrices $F$ and $V$ for the new infectious terms and the remaining transition terms, evaluated at the disease-free equilibrium are respectively given by

$F=\left[\begin{array}{cccc}0& {\beta}_{n1}& 0& 0\\ 0& 0& 0& 0\\ 0& \frac{{\Lambda}_{h}{\beta}_{n2}{\mu}_{n}}{{\Lambda}_{n}\left({\alpha}_{h}+{\mu}_{h}\right)}& 0& \frac{{\beta}_{h}{\mu}_{h}}{\left({\alpha}_{h}+{\mu}_{h}\right)}\\ 0& 0& 0& 0\end{array}\right]$ (23)

and

$V=\left[\begin{array}{cccc}\left({\mu}_{n}+{\nu}_{n}\right)& -{\beta}_{n1}& 0& 0\\ -{\nu}_{n}& \left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)& 0& 0\\ 0& 0& \left({\mu}_{h}+{\nu}_{h}\right)& 0\\ 0& 0& -{\nu}_{h}& \left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right)\end{array}\right]$ (24)

Therefore;

$F{V}^{-1}=\left[\begin{array}{cccc}\frac{{\nu}_{n}{\beta}_{n1}}{{y}_{n}}& \frac{{\beta}_{n1}}{\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)}& 0& 0\\ 0& 0& 0& 0\\ \frac{{\Lambda}_{h}{\beta}_{n2}{\nu}_{n}{\mu}_{n}}{{\Lambda}_{n}{y}_{n}\left({\alpha}_{h}+{\mu}_{h}\right)}& \frac{{\Lambda}_{h}{\beta}_{n2}{\mu}_{n}}{{\Lambda}_{n}\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)\left({\alpha}_{h}+{\mu}_{h}\right)}& \frac{{\beta}_{h}{\nu}_{h}{\mu}_{h}}{{y}_{h}\left({\alpha}_{h}+{\mu}_{h}\right)}& \frac{{\beta}_{h}{\mu}_{h}}{\left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right)\left({\alpha}_{h}+{\mu}_{h}\right)}\\ 0& 0& 0& 0\end{array}\right]$ (25)

where ${y}_{n}=\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)\left({\mu}_{n}+{\nu}_{n}\right)$ , ${y}_{h}=\left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right)\left({\mu}_{h}+{\nu}_{h}\right)$ .

Hence, the basic reproduction numbers of the model are given by:

${R}_{0}=\left\{{R}_{0,n},{R}_{0,h}\right\}$

where ${R}_{0,n}$ and ${R}_{0,h}$ are the monkeypox induced reproduction numbers for non-human primates and humans respectively and are given as:

${R}_{0,n}=\frac{{\nu}_{n}{\beta}_{n1}}{\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)\left({\mu}_{n}+{\nu}_{n}\right)}$ (26)

${R}_{0,h}=\frac{{\nu}_{h}{\beta}_{h}{\mu}_{h}}{\left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right)\left({\mu}_{h}+{\nu}_{h}\right)\left({\alpha}_{h}+{\mu}_{h}\right)}$ (27)

Theorem 1: The disease-free equilibrium is locally asymptotically stable if ${R}_{0}<1$ , and unstable if ${R}_{0}>1$ with ${R}_{0}=\mathrm{max}\left\{{R}_{0,n},{R}_{0,h}\right\}$ .

3.2.2. Local Stability of the Endemic Equilibrium (EE) Point

The local stability will be established using linearization method. Therefore, the Jacobian matrix $J$ of the model equations is given as:

$J=\left[\begin{array}{ccccccccc}-\left({\mu}_{n}+x\right)& 0& -m& 0& 0& 0& 0& 0& 0\\ x& -\left({\mu}_{n}+{\nu}_{n}\right)& m& 0& 0& 0& 0& 0& 0\\ 0& {\nu}_{n}& -{j}_{n}& 0& 0& 0& 0& 0& 0\\ 0& 0& {\rho}_{n}& -{\mu}_{n}& 0& 0& 0& 0& 0\\ 0& 0& -n& 0& -\left({\mu}_{h}+{\alpha}_{h}+w\right)& 0& 0& -z& 0\\ 0& 0& 0& 0& {\alpha}_{h}& -{\mu}_{h}& 0& 0& 0\\ 0& 0& n& 0& w& 0& -\left({\mu}_{h}+{\nu}_{h}\right)& 0& 0\\ 0& 0& 0& 0& 0& 0& {\nu}_{h}& -{j}_{h}& 0\\ 0& 0& 0& 0& 0& 0& 0& {\rho}_{h}& -{\mu}_{h}\end{array}\right]$ (28)

where $x=\frac{{\beta}_{n1}{I}_{n}^{\ast}}{{N}_{n}^{\ast}}$ , $m=\frac{{\beta}_{n1}{S}_{n}^{\ast}}{{N}_{n}^{\ast}}$ , $n=\frac{{\beta}_{n2}{S}_{h}^{\ast}}{{N}_{n}^{\ast}}$ , $w=\frac{{\beta}_{n2}{I}_{n}^{\ast}}{{N}_{n}^{\ast}}+\frac{{\beta}_{h}{I}_{h}^{\ast}}{{N}_{h}^{\ast}}$ , $z=\frac{{\beta}_{h}{S}_{h}^{\ast}}{{N}_{h}^{\ast}}$ ,

${j}_{h}=\left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right)$ and ${j}_{n}=\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)$ .

Next, we used elementary row-operations as used by [15] and [16] to row- reduce (28) to an upper triangular matrix and obtained the following eigenvalues:

${\psi}_{1}=-\left({\mu}_{n}+x\right)$ (29)

${\psi}_{2}=-\left({\mu}_{n}+x\right)\left({\mu}_{n}+{\nu}_{n}\right)$ (30)

${\psi}_{3}=-{r}_{n}$ (31)

${\psi}_{4}=-{r}_{n}{\mu}_{n}$ (32)

${\psi}_{5}=-{r}_{n}\left({\mu}_{h}+{\alpha}_{h}+w\right)$ (33)

${\psi}_{6}=-{\mu}_{h}{r}_{n}\left({\mu}_{h}+{\alpha}_{h}+w\right)$ (34)

${\psi}_{7}=-{r}_{n}^{2}\left({\mu}_{h}+{\alpha}_{h}+w\right)\left({\mu}_{h}+{\nu}_{h}\right)$ (35)

${\psi}_{8}=-{j}_{h}{r}_{n}^{2}\left({\mu}_{h}+{\alpha}_{h}+w\right)\left({\mu}_{h}+{\nu}_{h}\right)$ (36)

${\psi}_{9}=-{\mu}_{h}{j}_{h}{r}_{n}^{2}\left({\mu}_{h}+{\alpha}_{h}+w\right)\left({\mu}_{h}+{\nu}_{h}\right)$ (37)

where ${r}_{n}=\left[{j}_{n}\left({\mu}_{n}+x\right)\left({\mu}_{n}+{\nu}_{n}\right)-{\mu}_{n}{\nu}_{n}m\right]$ .

Therefore, since the real part of all the eigenvalues ${\psi}_{i}$ , for $i=1,2,\cdots ,9$ are negative, the endemic equilibrium is locally asymptotically stable from the following theorem:

Theorem 2: The endemic equilibrium is locally asymptotically stable if ${R}_{0}<1$ , and unstable if ${R}_{0}>1$ with ${R}_{0}=\mathrm{max}\left\{{R}_{0,n},{R}_{0,h}\right\}$ .

3.3. Global Stability Analysis of the Disease-Free Equilibrium Point

Theorem 3: The disease-free equilibrium is globally asymptotically stable if ${R}_{0}<1$ and unstable if ${R}_{0}>1$

Proof: By the comparison theorem, the rate of change of the variables representing the infectious classes in the model can be compared in the following inequality:

$\left[\begin{array}{c}{{E}^{\prime}}_{n}\\ {{I}^{\prime}}_{n}\\ {{E}^{\prime}}_{h}\\ {{I}^{\prime}}_{h}\end{array}\right]\le \left(F-V\right)\left[\begin{array}{c}{E}_{n}\\ {I}_{n}\\ {E}_{h}\\ {I}_{h}\end{array}\right]-{M}_{1}{\theta}_{1}{\left[\begin{array}{c}{E}_{n}\\ {I}_{n}\\ {E}_{h}\\ {I}_{h}\end{array}\right]}_{1}-{M}_{2}{\theta}_{2}\left[\begin{array}{c}{E}_{n}\\ {I}_{n}\\ {E}_{h}\\ {I}_{h}\end{array}\right]-{\theta}_{3}\left[\begin{array}{c}{E}_{n}\\ {I}_{n}\\ {E}_{h}\\ {I}_{h}\end{array}\right]$ (38)

where $F$ and $V$ are defined in (23) and (24) respectively, ${M}_{1}=1-\frac{{S}_{h}^{0}}{{N}_{h}^{0}}$ , ${M}_{2}=1-\frac{{V}_{h}^{0}}{{N}_{h}^{0}}$ , ${\theta}_{1},\text{\hspace{0.17em}}{\theta}_{2}$ and ${\theta}_{3}$ are nonnegative matrices. And since ${S}_{h}^{0}\le {N}_{h}^{0}$ ,

then ${V}_{h}^{0}\le {N}_{h}^{0}$ . Therefore, from (38) we get:

$\left[\begin{array}{c}{{E}^{\prime}}_{n}\\ {{I}^{\prime}}_{n}\\ {{E}^{\prime}}_{h}\\ {{I}^{\prime}}_{h}\end{array}\right]\le \left(F-V\right)\left[\begin{array}{c}{E}_{n}\\ {I}_{n}\\ {E}_{h}\\ {I}_{h}\end{array}\right]$ (39)

Therefore the matrix $\left(F-V\right)$ is obtained as:

$\left(F-V\right)=\left[\begin{array}{cccc}-\left({\mu}_{n}+{\nu}_{n}\right)& {\beta}_{n1}& 0& 0\\ {\nu}_{n}& -\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)& 0& 0\\ 0& \frac{{\Lambda}_{h}{\beta}_{n2}{\mu}_{n}}{{\Lambda}_{n}\left({\alpha}_{h}+{\mu}_{h}\right)}& -\left({\mu}_{h}+{\nu}_{h}\right)& \frac{{\beta}_{h}{\mu}_{h}}{\left({\alpha}_{h}+{\mu}_{h}\right)}\\ 0& 0& {\nu}_{h}& -\left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right)\end{array}\right]$ (40)

From the matrix in (40), let $\psi $ be an eigenvalue. Then, the characteristic equation $\left|\left(F-V\right)-\psi I\right|=0$ gives the following eigenvalues:

${\psi}_{10}=-\left[\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)\left({\mu}_{n}+{\nu}_{n}\right)-{\beta}_{n1}{\nu}_{n}\right]$ (41)

${\psi}_{11}=-\left({\mu}_{n}+{d}_{n}+{\rho}_{n}\right)$ (42)

${\psi}_{12}=-\left[\left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right)\left({\mu}_{h}+{\nu}_{h}\right)-\frac{{\beta}_{h}{\nu}_{n}{\mu}_{h}}{\left({\alpha}_{h}+{\mu}_{h}\right)}\right]$ (43)

${\psi}_{13}=-\left({\mu}_{h}+{d}_{h}+{\rho}_{h}\right)$ (44)

Therefore, all the four eigenvalues of the matrix in (40) have negative real part, showing that the matrix (40) is stable if ${R}_{0}<1$ . Consequently, using the model equations in (1)-(13), $\left({E}_{n},{I}_{n},{E}_{h},{I}_{h}\right)\Rightarrow \left(0,0,0,0\right)$ as $t\Rightarrow \infty $ . Thus by the comparison theorem as used in [17] , $\left({E}_{n},{I}_{n},{E}_{h},{I}_{h}\right)\Rightarrow \left(0,0,0,0\right)$ as $t\Rightarrow \infty $ . Evaluating the model system (3)-(11) at ${E}_{n}={I}_{n}={E}_{h}={I}_{h}=0$ gives

${S}_{h}^{0}=\frac{{\Lambda}_{h}}{\left({\alpha}_{h}+{\mu}_{h}\right)}$ , ${V}_{h}^{0}=\frac{{\Lambda}_{h}}{{\mu}_{h}}\frac{{\alpha}_{h}}{\left({\alpha}_{h}+{\mu}_{h}\right)}$ , ${S}_{n}^{0}=\frac{{\Lambda}_{n}}{{\mu}_{n}}$ and $\left({R}_{n},{R}_{h}\right)\Rightarrow \left(0,0\right)$ as

$t\Rightarrow \infty $ for ${R}_{0}<1$ . Hence, the disease-free equilibrium is globally asymptotically stable for ${R}_{0}<1$ .

4. Numerical Simulations and Sensitivity Analysis of Parameters

4.1. Numerical Simulations for the Model

In this section, numerical simulations for the model were carried out using the parameter values in Table 1. Some of these parameters were sourced from existing literatures where available, and assumed for the purpose of illustrations to fit the model analysis where otherwise. We used MATLAB R2012b encoded with ODE45 solver to simulate the model system using the parameters and an initial population of ${S}_{n}=250$ , ${E}_{n}=125$ , ${I}_{n}=75$ , ${R}_{n}=50$ , ${S}_{h}=8000$ , ${V}_{h}=5000$ , ${E}_{h}=3000$ , ${I}_{h}=2000$ and ${R}_{h}=2000$ .

4.2. Sensitivity Analysis of Parameters in the Model

Sensitivity indices allow us to measure the relative change in a variable when a parameter changes. The normalized forward sensitivity index of a variable to a parameter is the ratio of the relative change in the variable to the relative change in the parameter [18] . When the variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives from the following:

Definition: The normalized forward sensitivity index of a variable $\tau $ that depends, differentially, on a parameter $p$ , is defined as:

${\Upsilon}_{p}^{\tau}=\frac{\partial \tau}{\partial p}\times \frac{p}{\tau}$ (45)

We computed the sensitivity index of each parameter involved in ${R}_{0}$ using the parameter values in Table 1.

The indices with positive signs show that the value of ${R}_{0}$ increases when the corresponding parameters are increased and indices with negative signs indicates that, the value of ${R}_{0}$ decreases with increase in the corresponding parameters. This analysis is done to ascertain which parameters dominate the results of our analysis. Therefore, some parameters are deliberately excluded out of the sensitivity analysis due to their relative low importance in the actual disease transmission process. For example, the natural births, deaths in both humans and the non-human primates. The results of the analysis are presented in Table 2.

Therefore, it is clear from the Table 2 above, that ${R}_{0}$ will decrease with increase in the values of the control parameters ${\alpha}_{h}$ and ${\rho}_{h}$ since the sensitivity indices of these parameters are negative.

5. Results and Discussions

5.1. Results

The results of the analysis for the model were presented in Section 3 of this paper. The results of the numerical simulations for the model and sensitivity analysis of the model parameters using parameter values in Table 1 were presented in Figure 2 and Table 2 of Section 4 respectively. The computed basic reproduction numbers for the model using parameter values in Table 1 were

${R}_{0,n}=3.375\times {10}^{-3}$ (46)

${R}_{0,h}=9.1304\times {10}^{-6}$ (47)

Table 2. Numerical values of sensitivity indices for model parameters in ${R}_{0}$ and ${\lambda}_{h}^{\ast}$ .

(a) (b) (c) (d)

(e) (f) (g) (h)

(i)

Figure 2. Results of simulations with vaccination and treatment interventions, where ${R}_{0,n}=3.375\times {10}^{-3}$ and ${R}_{0,h}=9.1304\times {10}^{-6}$ .

Clearly, from (45) and (46), ${R}_{0,n}<1$ and ${R}_{0,h}<1$ , suggesting that the disease-free equilibrium is both locally and globally asymptotically stable while the endemic equilibrium of the model is locally asymptotically stable from our analysis.

The sensitivity indices of the model parameters in Table 2 evaluated using the parameter values in Table 1 suggest that, the indices with positive signs increases the value of ${R}_{0}$ when the corresponding parameters are increased and indices with negative signs decreases the value of ${R}_{0}$ with increase in the corresponding parameters.

5.2. Discussions

In this paper, we studied the dynamics of the transmission of the monkeypox virus infection under the combined vaccine and treatment interventions using the work of [9] as frame. An additional compartment representing the Latent or Exposed populations of the non-human primates and the humans was added to the existing work by [9] due to the identified fact in [5] , that monkeypox virus has incubation rates. As seen in the model diagram, the vaccine was administered on the susceptible human population with the assumption it confers permanent immunity against monkeypox virus infection at an initial rate of ${\alpha}_{h}=0.1$ . This low vaccination rate was used due to the identified increase in the prevalence of HIV/AIDS in its endemic environment as a result. However, the sensitivity analysis carried out in the model revealed that, if the control parameters of the treatment ${\rho}_{h}$ and vaccination ${\alpha}_{h}$ are increased, the basic reproduction numbers of the model which serves as the threshold for measuring new monkeypox virus infections among the two interacting populations will decrease.

The results of the numerical simulations carried out for the model using parameter values in Table 1 shows that, the infectious classes in both non-human primates and humans will be wiped out in the time considered by this study, whereby each of the infectious population becomes asymptotic to zero. This mean that, the disease-free equilibrium as seen in Figure 2, is asymptotically stable both locally and globally. Figure 2(e) and Figure 2(f) shows the susceptible human population decreasing exponentially while the vaccinated human population was growing exponentially up to equilibrium level before it started decreasing respectively. This can be explained as; due to the administration of the vaccine, the susceptible human population will continue to decrease resulting into most of the individuals in the class being vaccinated. Besides, the susceptible humans also suffer natural mortality. While increase in the vaccinated class can be explained from the continuous vaccination being carried out on the susceptible humans and the decrease in the population as seen in the Figure 2(f) was due to the fact, the compartment was only recharged by the low vaccination rate ${\alpha}_{h}$ and the class also suffers natural mortality. This fact suggests that, the vaccine rate can be increased for greener results and the vaccine should be re-administered whenever there is an outbreak of the virus in the future.

The treatment intervention as seen in Figure 2(i) caused the recovered class to grow exponentially up to equilibrium level, and which then started decreasing. This is due to the fact that the recovered human population is recharged by treating the infected humans. And this means that, when the infected human population approaches zero, the recovered class dies out exponentially, and besides, humans recover with permanent immunity and that recovered class also suffers natural mortality. As seen in Figure 2(d), the recovered non-humans class grows exponentially up to equilibrium and then dies out exponentially in the absence of an infected non-human primate. This population becomes asymptotic to zero with a smooth curve.

6. Conclusion

In this paper, we developed a mathematical model for the dynamics of transmission of the monkeypox virus infection with combined interventions of vaccination and treatment. We carried analysis on the developed model. The disease-free equilibrium was found to be both locally and globally asymptotically stable if ${R}_{0}<1$ and unstable if ${R}_{0}>1$ . Using parameter values obtained from existing literatures, we carried out numerical simulations and sensitivity analysis for the model and the parameters respectively. The simulations results revealed that, the disease will be eradicated from both humans and the non-human primates with the proposed interventions of the model in due time. Sensitivity analysis revealed that, the interventions offer an optimal control on the monkeypox virus infection in the human population with increase in the control parameter rates of vaccination and treatment.

References

[1] CDC (2003) Update: Multistate Outbreak of Monkeypox—Illinois, Indiana, Kansas, Missouri, Ohio and Wisconsin.

https://www.cdc.gov/mmwr/preview/mmwrhtml/mm5227a5.htm

[2] Jezek, Z., Szczeniowski, M., Paluku, K.M., Mutombo, M. and Grab, B. (1988) Human Monkeypox: Confusion with Chickenpox. Acta Tropica, 45, 297-307.

[3] Ladnyj, I.D., Ziegler, P. and Kima, E. (1972) A Human Infection Caused by Monkeypox Virus in Basankusu Territory, Democratic Republic of Congo. Bulletin of the World Health Organization, 46, 593-597.

[4] Von Magnus, P., Andersen, E.K., Petersen, K.B. and Birch-Anderson, A.A. (1959) A Pox-Like Disease in Cynomolgus Monkeys. Acta Pathol Microbiol Immunol Scand, 46, 156-176.

https://doi.org/10.1111/j.1699-0463.1959.tb00328.x

[5] CDC (2003) What You Should Know about Monkeypox.

https://www.cdc.gov/poxvirus/monkeypox/

[6] Sola, O. (2017) NCDC Confirms 12 Suspected Monkeypox Cases in Bayelsa. Vanguard Online Newspaper.

[7] The Eagle Online (2017) Breaking: FG Confirms Only Four Cases.

https://theeagleonline.com.ng/breaking-monkey-pox-fg-confirms-only-four-cases/

[8] World Health Organization (WHO) (1999) Technical Advisory Group on Human Monkeypox Report of a WHO Meeting. WHO, Geneva.

http://www.who.int/emc

[9] Bhunu, C.P. and Mushayabasa, S. (2011) Modeling the Transmission Dynamics of Pox-Like Infections. International Journal of Applied Mathematics, 41, 2.

[10] Niwas, L. (3003) Mathematical Modeling of Smallpox with Optimal Intervention Policy. Master’s Thesis, University of Central Florida, Orlando.

[11] Babak, P., Lauren, A.M., Danuta, M.S., Mel, K., David, M.P. and Robert, C.B. (2005) Modeling Control Strategies of Respiratory Pathogens. Emerging Infectious Diseases, 11, 1249-1256.

http://www.cdc.gov/eid

[12] Birkhof, G. and Rota, G.C. (1989) Ordinary Differential Equations. MIT Press, Boston.

[13] Van den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48.

https://doi.org/10.1016/S0025-5564(02)00108-6

[14] Diekmann, O., Heesterbeek, J.A. and Metz, J.A.J. (1990) On the Definition and the Computation of the Basic Reproductive Ratio, R0 in Models of Infectious Diseases in Heterogeneous Populations. Journal of Mathematical Biology, 28, 365-382.

https://doi.org/10.1007/BF00178324

[15] Usman, S., Adamu, I.I. and Aliyu, H.B. (2017) Mathematical Model for the Transmission Dynamics of Zika Virus Infection with Combined Vaccination and Treatment Interventions. Journal of Applied Mathematics and Physics, 5, 1964-1978.

[16] Usman, S., Adamu, I.I. and Dahiru, U. (2017) Stability Analysis of a Mathematical Model for the Transmission Dynamics of Zika Virus Infection. Journal of the Nigerian Association of Mathematical Physics, 40.

https://www.researchgate.net/publication/320617743

[17] Shaban, N. and Hawa, M. (2014) Modeling the Impact of Vaccination and Screening on the Dynamics of Human Papillomavirus Infection. International Journal of Mathematical Analysis, 8, 441-454.

http://dx.doi.org/10.12988/ijma.2014.312302

https://doi.org/10.12988/ijma.2014.312302

[18] Chitnis, N., Hyman, J.M. and Cushing, J.M. (2008) Determining Important Parameters in the Spread of Malaria through the Sensitivity Analysis of a Mathematical Model. Bulletin of Mathematical Biology, 70, 1272.

https://doi.org/10.1007/s11538-008-9299-0