Modeling the Transmission Dynamics of the Monkeypox Virus Infection with Treatment and Vaccination Interventions ()
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
(Sn) , Exposed/Latent
(En) , Infected
(In) and Recovered
(Rn) subpopulations. The non-human primates and/or some wild animals are recruited into the Susceptible
(Sn) class at a constant birth rate
Λn and become exposed to the monkeypox virus after getting into contact with an infected non-human primate at a rate
λn with
λn=βn1InNn (1)
where
β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
(In) at a rate
νn . The infected animals in
(In) are capable of either; infecting other animals when they come into contact, die due to the disease at a rate
dn or recover naturally with permanent immunity at a rate
ρn and move into
(Rn) . All non-human primates in the model experience natural mortality rate
μn .
The total human host population was also divided into Susceptible
(Sh) , Vaccinated
(Vh) , Exposed/Latent
(Eh) , Infected
(Ih) and Recovered
(Rh) human subpopulations. The Susceptible humans are recruited into
(Sh) through birth and migration at a constant rate
Λh . A susceptible individual is either vaccinated against the monkeypox virus at a rate
αh and move to
(Vh) 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
λh with
λh=βn2InNn+βhIhNh (2)
where
β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
β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
(Eh) proceeds to the infected class
(Ih) at a rate
νh . Individuals in
(Ih) either die due to the virus at a constant rate
dh or recover with permanent immunity after receiving treatment at a rate
ρh into
(Rh) . All individuals in the human subpopulations suffer a natural mortality at a constant rate
μ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′n=Λn−(μn+λn)Sn (3)
E′n=λnSn−(μn+νn)En (4)
I′n=νnEn−(μn+dn+ρn)In (5)
Parameter | Value | Source |
| 0.2 | Assumed |
| 0.029 | [7] |
| 0.1 | Assumed |
| 0.02 | [7] |
| 0.2 | Assumed |
| 0.1 | [7] |
| 0.3 | Assumed |
| 0.83 | [7] |
| 0.3 | Assumed |
| 0.095 | [4] |
| 0.1 | Assumed |
| 0.0027 | [7] |
| 0.00252 | [7] |
| 0.000063 | [7] |
翻译:
Table 1. Model parameter values.
R′n=ρnIn−μnRn (6)
S′h=Λh−(μh+λh+αh)Sh (7)
V′h=αhSh−μhVh (8)
E′h=λhSh−(μh+νh)Eh (9)
I′h=νhEh−(μh+dh+ρh)Ih (10)
R′n=ρhIh−μhRh (11)
Nn(t)=Sn+En+In+Rn (12)
Nh(t)=Sh+Vh+Eh+Ih+Rh (13)
Subject to the following nonnegative initial conditions:
Sn(0)≥0,En(0)≥0, In(0)≥0, Rn(0)≥0 (14)
Sh(0)≥0,Vh(0)≥0,Eh(0)≥0, Ih(0)≥0, Rh(0)≥0 (15)
Sn(0)+En(0)+In(0)+Rn(0)≤Nn(0) andSh(0)+Vh(0)+Eh(0)+ Ih(0)+ Rh(0)≤Nh(0) (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
Ω . Thus the feasible region
Ω={(Sn,En,In,Rn)∈ℝ4+ : Nn≤Λnμn(Sh,Vh,Eh,Ih,Rh)∈ℝ5+ : Nh≤Λhμh (17)
is considered. Therefore, after differentiation of (12) and (13), and proper substitutions, we have:
dNn(t)dt=Λn−μnNn−dnIn≤Λn−μnNn (18)
and;
dNh(t)dt=Λh−μhNh−dhIh≤Λh−μhNh (19)
Applying [12] on the differential inequalities in (18) and (19), we obtained:
{Nn(t)≤Nn(0)e−μnt+Λnμn(1−e−μnt)Nh(t)≤Nh(0)e−μht+Λhμh(1−e−μht) (20)
where
Nn(0) and
Nh(0) are the initial populations of the non-human
primates and the humans respectively. Therefore,
0≤Nn≤Λnμn and
0≤Nh≤Λhμh as
t→∞ . This implies that,
Λnμn and
Λhμh are upper bounds for
Nn(t) and
Nh(t) respectively, as long as
Nn(0)≤Λnμn and
Nh(0)≤Λhμh . Hence, the
feasible solution of the model equations in (3)-(13) enters the region
Ω which is a positively invariant set. Thus, the system is mathematically and epidemiologically well-posed. Therefore, for an initial starting point
x∈Ω , the trajectory lies in
Ω , and so it is sufficient to restrict our analysis on
Ω . Clearly, under the dynamics described by the model equations, the closed set
Ω is hence a positively invariant set.
3.1. Model Equilibrium Points
Using standard approaches, the model disease-free
ε0 and endemic
ε∗ (which existed only when
R0>1 ) equilibrium points are established as follows:
ε0=(Sh,Vh,Eh,Ih,Rh,Sn,En,In,Rn)=(Λhαh+μh,Λhμhαhαh+μh,0,0,0,Λnμn,0,0,0) (21)
ε∗=(S∗h,V∗h,E∗h,I∗h,R∗h,S∗n,E∗n,I∗n,R∗n) (22)
where:
S∗h=Λh(μh+αh+λ∗h) ,
V∗h=Λhμhαh(μh+αh+λ∗h) ,
E∗h=λ∗hΛh(μh+νh)(μh+αh+λ∗h) ,
I∗h=νhλ∗hΛh(μh+dh+ρh)(μh+νh)(μh+αh+λ∗h) ,
R∗h=Λhμhρhνhλ∗h(μh+dh+ρh)(μh+νh)(μh+αh+λ∗h) ,
S∗n=Λn(μn+λ∗n) ,
E∗n=λ∗nΛn(μn+νn)(μn+λ∗n) ,
I∗n=νnλ∗nΛn(μn+νn)(μn+λ∗n)(μn+dn+ρn) ,
R∗n=Λnμnρnνnλ∗n(μn+νn)(μn+λ∗n)(μn+dn+ρn) with
λ∗n=βn1I∗nN∗n ,
λ∗h=βn2I∗nN∗n+βhI∗hN∗h ,
N∗n=Λn−dnI∗nμn and
N∗h=Λh−dhI∗hμ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
|FV−1−ψI|=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=[0βn10000000Λhβn2μnΛn(αh+μh)0βhμh(αh+μh)0000] (23)
and
V=[(μn+νn)−βn100−νn(μn+dn+ρn)0000(μh+νh)000−νh(μh+dh+ρh)] (24)
Therefore;
FV−1=[νnβn1ynβn1(μn+dn+ρn)000000Λhβn2νnμnΛnyn(αh+μh)Λhβn2μnΛn(μn+dn+ρn)(αh+μh)βhνhμhyh(αh+μh)βhμh(μh+dh+ρh)(αh+μh)0000] (25)
where
yn=(μn+dn+ρn)(μn+νn) ,
yh=(μh+dh+ρh)(μh+νh) .
Hence, the basic reproduction numbers of the model are given by:
R0={R0,n,R0,h}
where
R0,n and
R0,h are the monkeypox induced reproduction numbers for non-human primates and humans respectively and are given as:
R0,n=νnβn1(μn+dn+ρn)(μn+νn) (26)
R0,h=νhβhμh(μh+dh+ρh)(μh+νh)(αh+μh) (27)
Theorem 1: The disease-free equilibrium is locally asymptotically stable if
R0<1 , and unstable if
R0>1 with
R0=max{R0,n,R0,h} .
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=[−(μn+x)0−m000000x−(μn+νn)m0000000νn−jn00000000ρn−μn0000000−n0−(μh+αh+w)00−z00000αh−μh00000n0w0−(μh+νh)00000000νh−jh00000000ρh−μh] (28)
where
x=βn1I∗nN∗n ,
m=βn1S∗nN∗n ,
n=βn2S∗hN∗n ,
w=βn2I∗nN∗n+βhI∗hN∗h ,
z=βhS∗hN∗h ,
jh=(μh+dh+ρh) and
jn=(μn+dn+ρn) .
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:
ψ1=−(μn+x) (29)
ψ2=−(μn+x)(μn+νn) (30)
ψ3=−rn (31)
ψ4=−rnμn (32)
ψ5=−rn(μh+αh+w) (33)
ψ6=−μhrn(μh+αh+w) (34)
ψ7=−r2n(μh+αh+w)(μh+νh) (35)
ψ8=−jhr2n(μh+αh+w)(μh+νh) (36)
ψ9=−μhjhr2n(μh+αh+w)(μh+νh) (37)
where
rn=[jn(μn+x)(μn+νn)−μnνnm] .
Therefore, since the real part of all the eigenvalues
ψi , for
i=1,2,⋯,9 are negative, the endemic equilibrium is locally asymptotically stable from the following theorem:
Theorem 2: The endemic equilibrium is locally asymptotically stable if
R0<1 , and unstable if
R0>1 with
R0=max{R0,n,R0,h} .
3.3. Global Stability Analysis of the Disease-Free Equilibrium Point
Theorem 3: The disease-free equilibrium is globally asymptotically stable if
R0<1 and unstable if
R0>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:
[E′nI′nE′hI′h]≤(F−V)[EnInEhIh]−M1θ1[EnInEhIh]1−M2θ2[EnInEhIh]−θ3[EnInEhIh] (38)
where
F and
V are defined in (23) and (24) respectively,
M1=1−S0hN0h ,
M2=1−V0hN0h ,
θ1, θ2 and
θ3 are nonnegative matrices. And since
S0h≤N0h ,
then
V0h≤N0h . Therefore, from (38) we get:
[E′nI′nE′hI′h]≤(F−V)[EnInEhIh] (39)
Therefore the matrix
(F−V) is obtained as:
(F−V)=[−(μn+νn)βn100νn−(μn+dn+ρn)000Λhβn2μnΛn(αh+μh)−(μh+νh)βhμh(αh+μh)00νh−(μh+dh+ρh)] (40)
From the matrix in (40), let
ψ be an eigenvalue. Then, the characteristic equation
|(F−V)−ψI|=0 gives the following eigenvalues:
ψ10=−[(μn+dn+ρn)(μn+νn)−βn1νn] (41)
ψ11=−(μn+dn+ρn) (42)
ψ12=−[(μh+dh+ρh)(μh+νh)−βhνnμh(αh+μh)] (43)
ψ13=−(μh+dh+ρh) (44)
Therefore, all the four eigenvalues of the matrix in (40) have negative real part, showing that the matrix (40) is stable if
R0<1 . Consequently, using the model equations in (1)-(13),
(En,In,Eh,Ih)⇒(0,0,0,0) as
t⇒∞ . Thus by the comparison theorem as used in [17] ,
(En,In,Eh,Ih)⇒(0,0,0,0) as
t⇒∞ . Evaluating the model system (3)-(11) at
En=In=Eh=Ih=0 gives
S0h=Λh(αh+μh) ,
V0h=Λhμhαh(αh+μh) ,
S0n=Λnμn and
(Rn,Rh)⇒(0,0) as
t⇒∞ for
R0<1 . Hence, the disease-free equilibrium is globally asymptotically stable for
R0<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
Sn=250 ,
En=125 ,
In=75 ,
Rn=50 ,
Sh=8000 ,
Vh=5000 ,
Eh=3000 ,
Ih=2000 and
Rh=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
τ that depends, differentially, on a parameter
p , is defined as:
ϒτp=∂τ∂p×pτ (45)
We computed the sensitivity index of each parameter involved in
R0 using the parameter values in Table 1.
The indices with positive signs show that the value of
R0 increases when the corresponding parameters are increased and indices with negative signs indicates that, the value of
R0 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
R0 will decrease with increase in the values of the control parameters
αh and
ρ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
R0,n=3.375×10−3 (46)
R0,h=9.1304×10−6 (47)
Parameter Symbol | Sensitivity Index |
| −0.133 |
| −0.089 |
| −0.15 |
| −0.105 |
| +0.25 |
| +0.17 |
| −0.83 |
| +1.00 |
| +0.168 |
| +1.00 |
翻译:
Table 2. Numerical values of sensitivity indices for model parameters in
R0 and
λ∗h .
Clearly, from (45) and (46),
R0,n<1 and
R0,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
R0 when the corresponding parameters are increased and indices with negative signs decreases the value of
R0 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
α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
ρh and vaccination
α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
α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
R0<1 and unstable if
R0>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.