Mathematical Model for the Dynamics of Neisseria Gonorrhea Disease with Natural Immunity and Treatment Effects

Neisseria gonorrhea infection; a sexually transmitted disease, is caused primarily by a type of germ; a bacteria called neisseria gonorrhea. The infection is a major public health challenge today due to the high incidence of infections accompanied by a dwindling number of treatment options especially in developing and underdeveloped countries. In this paper, we developed a mathematical model for the transmission dynamics of neisseria gonorrhea infection and studied the effect of natural immunity and treatment as the only available control interventions on the spread of the disease in a population. We computed the model disease-free equilibrium and analyzed its local and global stability in a well-defined positively invariant and attracting set  using the next-generation matrix plus linearization method and the comparison theorem respectively. The disease-free equilibrium was proved to be both locally and globally asymptotically stable if 0 1 R  and unstable if 0 1 R  . We conducted sensitivity analysis of parameters in the basic reproduction number 0 R using the normalized forward sensitivity index method. Results of the analysis revealed that 0 R decreases with increase in treatment and natural immunity rates. The results of the numerical simulations carried out using MATLAB R2012B showed that there is increase in new infections due to increased contact with infected individuals in the susceptible population and that, with increased treatment rate and controlled death due to the disease in the population, neisseria gonorrhea infection would be wiped out within 300 days of the treatment intervention.


Introduction
Gonorrhea is a major public health challenge today, due to the high incidence of infections accompanied by a dwindling number of treatment options (Schiffert Health Center, 2011).It is a sexually transmitted disease caused by a type of germ, bacteria called Neisseria gonorrhea.It is passed from one person to another during vaginal, anal and oral sex (Amalia, G. & Loukas, Z., 2015).Gonorrhea is one of the sexually lethal transmitted diseases (STD) due to the number of complications that it causes in the infected persons (Sacrifice, N. K., et al., 2016;Riley, S., 2010).The disease has an incubation period of 10 days in women and 2-5 days in men (Balmelli, C. & Guntlard, H. F., 2003), and it was estimated that, approximately 106 million new cases of gonorrhea among adults globally was estimated (Benedek, T. G., 2005).The bacteria can be found in the throat, vagina, Urethra, and anus (Mushayabasa, S., et al., 2010).The history of the discovery of the neisseria Gonorrhea infection can be traced to Edinburg in the year 1792, where the surgeon Benjamin Bell clearly differentiated it from syphilis infection (Benedek, T. G., 2005).Gonorrhea causes a lot of complications in the infected persons, and a prolong infection can lead to severe eye infections, infertility in both men and woman, ectopic pregnancy, spontaneous abortion, still births and eventually death if untreated (Riley, S., 2010).The bacteria targets the cells of the mucous membranes, including the surfaces of the Urethra, vagina cervix, anus, rectum, the lining of the eyelid and the throat (Schiffert Health Center, 2011).Behavioral change is very important in the control of the spread of diseases and as such, can be involved in the interpretation of the diseases outbreak data and estimation of decreases transmission rates (Prabhakararao, G., 2013).
The epidemiology of infectious diseases has moved beyond identifying etiological agents and risk factors to a more detailed understanding of the mechanisms for controlling the distribution of infections (Garnette, G. P., 2002).Mathematical models provide an explicit framework to study the transmission dynamics of infectious diseases (Garnette, G. P., 2002).In 2013, Prabhakararao developed a mathematical model for the transmission dynamics of gonorrhea infection in Anantapur District, Andhra-Pradesh-India.In developing the model, they partitioned the population into susceptible, infectious and removed compartments and assumed recovery from infection confers permanent immunity.
The result of their study showed that the epidemic does not die out, but ultimately approach a steady state with reference to its severity among the population of Anantapur District.A more robust study was carried out by Ana, L. & James, A. Y. in 1976, where they took into consideration the fact that; the spread of gonorrhea in a population is nonuniform, and by partitioning the population into n groups they developed a deterministic model for Gonorrhea in a non-homogenous population.They carried out a study of the asymptotic stability properties of the model, to show that the disease will die out either for all positive initial disease levels or for none, depending on the contact rates and the lengths of infectious periods.Furthermore, they showed that if the initial number of infectives in at least one group is not zero, then the disease will remain endemic and the number of infectives in each group will approach a constant positive level.They concluded that; this behavior of gonorrhea epidemics is different from most other infectious disease, which usually have occasional major outbreaks followed by low disease levels, and that this different behavior of gonorrhea is attributable to the lack of developed immunity and the basic non seasonality of the disease.
In their work, (Amalia, G. & Loukas, Z., 2015) they developed and investigated a deterministic epidemic model for the spread of gonorrhea.They showed that the discrete time dynamical system exhibits far more complex dynamics than its continuous analogues.They also carried out stability analysis of their model and established that there are phenomena of Fold Flip bifurcations.The result of their study revealed that male Latex condom use stabilizes the chaotic vibrations of the system to a point where the number of infected individuals remains stable and is significantly small or zero, leading to the control of the disease.Some studies looked at the impact of treatment on the dynamics of neisseria gonorrhea infection.Amongst them is the work by (Sacrifice, N. K., et al., 2016), where they developed a mathematical model for the dynamics of gonorrhea with treatment effect and studied the model dynamics to understand the epidemic phenomenon and recommended strategies for its control.The result of their study showed that an increase in treatment rate has significant effect on the infectives.However, in developing their model equations, they III) Treated individuals move to exposed compartment.In reality, treated individuals acquire susceptible status to move to susceptible compartment or die naturally, its only after interacting with an infectious individual that a treated individual becomes exposed since recovery from gonorrhea does not confer permanent immunity.Again, incorporating this term in the model equations is not correct In this work, we reviewed the existing model by Sacrifice, N. K., et al., 2016, and developed a more realistic model that addresses the issues raised, by incorporating natural immunity effects into the dynamics of neisseria gonorrhea infection.

Model Assumptions
The model in this paper was developed under the following specific assumptions:  The population is finite with a constant recruitment rate.
 Members of the population interact freely not subject to any quarantine policy  The population is recharged by birth only.
 No emigration or immigration was considered in the model, hence the population is closed.
 Condom use and other protective measures were not considered in this model.
 Recovered individual becomes susceptible since recovery from gonorrhea does not confer permanent immunity.
 All members in the population suffer natural mortality at a uniform rate.
 Individuals in the subpopulations exhibit constant natural immunity response against the disease.
 Model parameters are strictly nonnegative and will be sourced from existing literatures or will be hypothetically estimated otherwise.

Definitions of State Variables and Parameters of the Model
The variables and parameters of the model are presented in Table 1.
Subject to the following nonnegative initial conditions:

Results
In this section, we begin the model analysis by showing that all feasible solutions of the model are uniformly bounded in a proper subset of  .The feasible region Therefore, solving the differential inequality in   .Hence, the feasible solution of the model system enters the region  which is a positively invariant set.Thus, the system is mathematically well-posed and epidemiologically meaningful.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.

Stability Analysis of the Disease-Free Equilibrium (DFE)
The disease-free equilibrium of the model was obtained by solving the model equations simultaneously with   0.

It was established as follows:
    0 , , , , , 0, 0, 0, 0 , 0, 0, 0, 0 It is important to notice that, in the absence of the disease (neisseria gonorrhea infection), the susceptible population tends to represents the entire host population, and is bounded by only birth and natural mortality rates of the population, i. e.,

Local Stability Analysis of the Disease-Free Equilibrium (DFE)
The local stability of the disease-free equilibrium follows immediately from the following theorem: Theorem1: The disease-free equilibrium is locally asymptotically stable if 0 1 R  and unstable if 0 1 R  .

Proof:
The basic reproduction number of the model was computed using the next-generation matrix as defined in Van den Driessche & Watmough (2002) and Diekmann, Heesterbeek & Metz (1990).It is defined to be largest eigenvalue or spectral radius of the characteristic equation 13 which is the largest eigenvalue of the characteristic polynomial equation.The proof of Theorem1 can be established by constructing a Jacobian matrix for the model system in order to analyze the local stability of the DFE.The Jacobian matrix evaluated at DFE wss obtained as The reduced Jacobian matrix in   16 gives the following eigenvalues Therefore, obviously 1,2 0   .But, 3,4,5 0   if 0 1 R  .Hence, the disease-free equilibrium (DFE) is locally asymptotically stable if 0 1 R  , which complete the proof of Theorem1.

Global Stability Analysis of the Disease-Free Equilibrium (DFE)
Theorem2: The disease-free equilibrium is globally asymptotically stable if 0 1 R  and unstable if 0 1 R  .
Proof: The proof of Theorem2 will be established using the comparison theorem as used Usman &Adamu 2017 andShaban et al. 2014.Thus, by the comparison theorem, the rate of change of the variables representing the infectious compartments in the model can be compared in the following inequality:  , we obtained the following eigenvalues: Therefore, all the eigenvalues of the matrix in   19 have negative real parts showing that the matrix in   , and     , 0,0 TR as t  for 0 1 R  .Hence, the disease-free equilibrium is globally asymptotically stable for 0 1 R  which complete the proof of Theorem2.

Numerical Simulations
In this section, numerical simulations for the model were carried out using the baseline parameter values in Appendix B, Table 2.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 parameter values and an assumed initial subpopulation of the human host.The results of the numerical simulations are given in Appendix A, Figure 2-10.In the simulations, we varied the values of the parameters 1 ,  and  accordingly to get the graphs.In each case, the value of the basic reproduction number 0 R was computed and displayed along with graphs.

R
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 as defined by Usman et al. 2017 is the ratio of the relative change in the variable to the relative change in the parameter.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, differentiably, on a parameter p , is defined as: We computed the sensitivity index of each parameter in 0 The indices with positive signs show that, the value of 0 R increases when the corresponding parameters are increased, and indices with negative signs indicates that, the value of 0 R decreases with increase in the corresponding parameters.
The results of the analyses are presented in Table 3, Appendix B.

Discussion
The results of the numerical simulations given in Figure A1, A2, A3, A4 and A5 of Appendix A shows that, if the parameter  is controlled at 2.0   , the disease will be wiped out of the population with the value of 0 R reaching as low as 0 0.10526 R  . But, when the contact rate is as high as 2.4

 
, even with treatment and natural immunity, the infectious classes reduces a little, then grows up to equilibrium point with asymptotic stability and does not die out over time.The results also revealed there is increase in new infections with increased contact rate  .Figure A6, A7, A8 and A9 shows the results of simulations with varying rates of; infection due to exposure, treatment and disease-induced death where  was kept at 2.0   .The basic reproduction numbers for the cases were also computed with the graphs.In each simulation, the four infectious classes die out exponentially and become asymptotic to zero within the first 300 days of the study with treatment intervention and the natural immunity rates kept constant.The susceptible class grows exponentially up to equilibrium level and achieved asymptotic stability in each case with a smooth curve.The basic reproduction number of the simulation result presented in Figure A8 and A9 coincide with that of the result in Figure 3, i.e. 0 0.10526 R  which is the least of all results.This is achieved by increased treatment rate and controlled death due to the neisseria gonorrhea infection.The sensitivity analysis of parameters in 0 R revealed that, increase in treatment rate and boosting the natural immunity of individuals will help in reducing the appearance of new cases of neisseria gonorrhea infection in the population.
In this paper, we reviewed the model by N. K., et al., 2016 and reformulated the model equations based on the three observations raised in section 1.0.We also incorporated natural immunity effects into the dynamics of the neisseria gonorrhea disease.This ushered in, an improved mathematical model for the transmission dynamics of the disease under the combined effects of treatment and natural immunity.We established the model disease-free equilibrium DFE and basic reproduction number 0 R .The DFE was proved to be both locally and globally asymptotically stable if 0 1 R  and unstable if 0 1 R  .We carried out sensitivity analysis of parameters in the 0 R and computed their respective sensitivity indices.The sensitivity analysis result is in Table B2 of Appendix B. The result of the analysis revealed that, increase in treatment rate and boosting the natural immunity of individuals will help in reducing the appearance of new cases of neisseria gonorrhea infection in the population.Furthermore, qualitative analysis carried out on the model with combined effects of natural immunity and treatment using the parameter values in Table B1 of Appendix B, showed that the disease would be wiped out of the population within the first 300 days of the treatment intervention and controlled death due to the disease.

Figure 1 .
Figure 1.Schematic flow diagram of the model the notations inVan den Driessche & Watmough (2002) for the model system in proportions, 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 used elementary row operations as used by Usman, Adamu & Aliyu 2017 and Usman & Adamu 2017 to reduce the Jacobian matrix to an upper triangular matrix of the form:

Table 1 .
Definitions of variables and parameters Et Number of susceptible individuals at time t Number of exposed individuals at time t Number of infected individuals at time t Number of treated individuals at time t Number of recovered individuals at time t Total population at time t It after the disease incubation period at a rate .An individual in   Tt at a constant rate 1  or die as a result of the neisseria gonorrhea infection at a rate  .A treated individual in   Tt recover with temporary immunity conferred by treatment and move to   Rt moves back to the susceptible compartment   St at a rate  after waning their temporary immunity.All individuals in the five subpopulations suffer natural mortality at a uniform constant rate  .The cycle continues in this manner.All model parameters in this paper are strictly nonnegative.The schematic flow diagram of the model described in this section is given in figure 1.
Nt at time t is divided into five subpopulations, viz; Susceptible individuals   St, Exposed individuals   Et, Infected individuals   Rt .Individuals are recruited into the susceptible population   St at a constant birth rate  .A susceptible individual in   St comes into contact with neisseria gonorrhea infection through an infectious individual in   It at a rate  and move to   Et.An individual in   rate  .All individuals in  

Table B2 .
Numerical Values of Sensitivity Indices for Model Parameters in 0