A co-infection model of dengue and leptospirosis diseases

In this paper an SIR deterministic mathematical model for co-infection of dengue and leptospirosis is proposed. We use a compartment model by using ordinary differential equations. The positivity of future solution of the model, the invariant region, and the stability of disease-free equilibrium point as well as endemic equilibrium point are studied. To study the stability of the equilibria, a basic reproduction number is obtained by using the next generation matrix. The robustness of the model is also investigated. To identify the effect of each parameter on the expansion or control of the diseases, sensitivity analysis is performed. The effects of treating dengue infected only, leptospirosis infected only, and co-infected individuals have been identified by using the numerical simulation. Therefore, increasing the rate of recovery and decreasing the contact rate of dengue, leptospirosis, and their co-infection have a great influence in fighting dengue, leptospirosis, and their co-infection in the community.


Introduction
Dengue is the world's fastest spreading mosquito transmitted infectious disease. An estimated 50 million dengue infections occur annually, and nearly 2.5 billion people live in regions with potential risk of dengue transmission [1][2][3]. Roughly half of the world's population probably lives in regions ecologically suitable for dengue transmission [3]. Dengue is transmitted to humans by mosquitoes of Aedes species, which thrive around the globe in tropical and subtropical urban centers [1,2]. The virus is transmitted by the infected female mosquito called Aedes aegypti [3]. There are four serotypes, DENV1-DENV4 [2,3]. On the other hand, a zoonotic bacterial disease leptospirosis is one of the most significant neglected tropical bacterial diseases and one of the leading zoonotic causes of morbidity in the world [1,4]. The most important reservoirs are small mammals, with large herbivores as additional important sources of infection. Pathogenic species of Leptospira were isolated from hundreds of mammalian species including bats and pinnipeds [5]. Leptospirosis is a global occurrence particularly in tropical and subtropical countries [5,6]. It mainly affects vulnerable populations, with an estimated global annual incidence of 1.03 million people and 58.900 deaths [4,6]. It is transmitted by contact with infected animal prod-ucts and with water, wet surfaces contaminated with animal urine [6]. Leptospirosis and dengue co-infections (LDCI) have been reported in many countries [7][8][9][10]. Co-infections occur due to simultaneous transmission of both these infections during the rainy season, and both of these infections are characterized by fever with myalgia.
Mathematical modeling has a great role in describing the dynamics of infectious diseases in a community. Several scholars have developed different models for dengue, leptospirosis, and their co-infection with other diseases to study their dynamics. For example, many mathematical models have been developed to understand dengue [11][12][13][14][15][16] and leptospirosis [17][18][19][20]. Also some co-infection models for dengue and chikungunya [21,22] and for dengue-Zika [23,24] have been proposed and analyzed. However, to the best of our knowledge, no one has investigated co-infection of dengue and leptospirosis. Therefore, in this paper we are interested in filling this gap.
The paper is organized as follows: Sect. 2 is devoted to the description and formulation of the model. In Sect. 3, the analysis of dengue only model, leptospirosis only model, and the co-infection model is performed. In Sect. 4 numerical simulations to give a better interpretation of the analytical results are reported. Our discussion and conclusion are presented in Sect. 5.

Model description and formulation
In this paper, we consider a deterministic mathematical model for co-infection of dengue and leptospirosis in a human population. The total population is divided into seven subclasses, which are susceptible population (S), dengue infectious (I d ), leptospirosis infectious (I l ), dengue and leptospirosis co-infectious (I dl ), dengue recovered (R d ), leptospirosis recovered (R l ), and dengue leptospirosis co-infectious recovered (R dl ). Each vector population of both diseases is also divided into two susceptible vectors of dengue (S vd ), infected vector of dengue (I vd ), susceptible vector of leptospirosis (S vl ), and infected vector of leptospirosis (I vl ).
Susceptible are recruited with the rate of through birth or immigration, and their number increases from individuals that come from sub-classes of dengue recovered, leptospirosis recovered, and co-infectious recovered by losing their temporary immunity at a rate of σ 1 , σ 2 , and σ 3 , respectively.
The entire susceptible population join I d compartment, when individuals can get dengue with contact rate of β 1 from a dengue only infected vectors. In a similar way, individuals can get leptospirosis by the contact rate of α 1 from a leptospirosis only infected or coinfected person with force of infection of leptospirosis α 1 (I l + I dl ) or by the contact rate α 2 from the infected vector of leptospirosis (I vl ) and join I l compartment.
Dengue only infected individuals also get an additional leptospirosis infection by the contact rate α 1 from a leptospirosis only infected and co-infected person with force of infection of leptospirosis α 1 (I l + I dl ) or by the contact rate α 2 from the infected vector of leptospirosis (I vl ) and join co-infected compartment (I dl ). On the other hand, leptospirosis only infected individuals get an additional dengue infection by the contact rate β 1 from the infected vector of dengue (I vd ) and join co-infected compartment (I dl ).
Dengue only infected individuals recover with the rate of γ 1 and join dengue only recovered compartment (R d ). In a similar way, leptospirosis only infected individuals also recover with the rate of γ 2 and join leptospirosis only recovered compartment (R l ). The co-infected compartment also recovers with ε rate, but those individuals either recover Rate of recovery from dengue disease ε Rate of recovery from dengue-leptospirosis co-infection diseases Mortality rate related to leptospirosis vector due to the disease only from dengue disease and join dengue only recovered compartment with probability of εη, or recover only from leptospirosis only disease and join leptospirosis only recovered compartment with probability of εk(1η), or recover from both diseases and join co-infected recover compartment with probability of g = ε(1k) (1η). In all compartments μ is a natural death rate of human population. Moreover, δ 1 is a dengue only caused death rate, and δ 2 is a leptospirosis caused death rate. and are the recruitment rates for vector population of dengue and leptospirosis, respectively. The leptospirosis and dengue vector populations have the natural death rate μ vl and μ vd , respectively. The description of parameters of the model are found in Table 1.
With regards to the above considerations, we have the compartmental flow diagram shown in Fig. 1. From the flow chart, the model will be governed by the following system of differential equations:

Qualitative analysis
In this section, for simplification of the work, we split the full model into sub-models, which are dengue and leptospirosis only models. The qualitative behavior of the submodel is studied first and the full model follows.

Dengue fever only model
To get this model from equation (1), we set I l = I dl = R l = R dl = 0, σ 2 = σ 3 = 0, I l + I dl = 0, and then we get (2)

Invariant region
To get an invariant region, which shows as the solution is bounded, the total human population of the model is N h = S + I d + R d . Differentiating N h both sides and substituting respective expressions of dS dt , dI d dt , and dR d dt from equation (2), we get If we are not considering death from dengue fever, then equation (3) becomes Solving equation (4), we get Similarly, for the vector population N vd = S vl + I vl , we have Therefore, all the solution set of (2) is bounded in D = D hd × D vd .

Positivity of solutions
which can be taken as after evaluating equation (6), we obtain Similarly, we obtain Therefore, all the solution sets are nonnegative for t ≥ 0.

Basic reproduction number
The basic reproduction number is an average number of secondary cases of infection when a single infectious individual is introduced into the total susceptible population. To obtain (R 0d ), the next generation matrix method that was formulated by [25] is used, and it is given by

Local stability of DFEP
Proof The Jacobian matrix of the dengue only model at DFEP is given by The eigenvalues of the matrix in (8), -μ, -(σ 1 + μ) and -μ vd are clearly negative. The other two eigenvalues can be obtained from the quadratic equation given by where Here, it follows from Routh-Hurwitz criteria ψ 1 > 0 and ψ 2 > 0. Thus, all the eigenvalues of the matrix in equation (9) have negative real parts if R 0d < 1. Hence DFEP is locally asymptotically stable.

Existence of endemic equilibrium point (EEP)
The endemic equilibrium point is denoted by E * = (S * , I * d , R * d , S * vd , I * vd ) and it occurs when the disease persists in the community.

Lemma 3.3 For R 0d > 1 a unique EEP E * exists and no EEP otherwise.
Proof To obtain it, we equate all the model equations in equation (2) to zero. Then we obtain where

Existence of backward bifurcation of dengue only model
We investigated the nature of the bifurcation by using the method introduced in [26], which is based on center manifold theory [26]. In order to apply theory, the following rearrangement and modification of variables are done on the submodel in equation (2). We let S = z 1 , I d = z 2 , R d = z 3 , S vd = z 4 , and I vd = z 5 . In addition, using vector notation z = (z 1 , given in the following, we chose β 1 as a bifurcation parameter and solve R 0d = 1, which leads to where The Jacobian matrix evaluated at DFEP is given by (z 1 = μ , z 2 = 0; z 3 = 0, z 4 = μ vd , z 5 = 0): The right eigenvector u = (u 1 , u 2 , u 3 , u 4 , u 5 ) T associated with the simple zero eigenvalue can be obtained as follows: Similarly, we obtain the left eigenvector 5 ) associated with the simple zero eigenvalue and given by Since the first, third, and fourth components of v are zero, we do not need the derivatives of f 1 , f 3 , and f 4 . From the derivatives of f 2 and f 5 , the only ones that are nonzero are the following: and all the other partial derivatives are zero. The direction of the bifurcation at R 0d = 1 is determined by the signs of the bifurcation coefficients a and b, obtained from the above partial derivatives, given respectively by By the fact that coefficient b is positive, it shows that the sign of coefficient a determines the local dynamics around DFEP for β 1 . It is clearly seen that system (11) undergoes backward bifurcation, and therefore we propose the following theorem. (1) exhibits backward bifurcation at R 0d = 1 and the disease-free equilibrium may not be globally stable.

Leptospirosis only model
To get this model from equation (1), we set I d = I dl = R d = R dl = 0, σ 1 = σ 3 = 0, and then we get

Invariant region
To get an invariant region, the total population of the model is N hl = S + I l + R l . Then the derivative of N hl with respect to time gives If we are not considering death from leptospirosis, then equation (15) becomes After solving equation (16), we get Similarly, for the vector population N vl = S vl + I vl , we have Therefore, all the solution set of the model in equation (14) is bounded in D 2 = D l × D vl . Proof To show solutions of the model are positive, first we take dS dt from equation (14):

Local stability of DFEP
Theorem 3.6 DFEP is locally asymptotically stable if R 0l < 1 and unstable if R 0l > 1.
Proof The Jacobian matrix at DFEP is given by The eigenvalues -μ, -(σ 2 + μ) and -μ vl are clearly negative. The other two eigenvalues can be obtained from the quadratic equation given by We applied Routh-Hurwitz criteria, and by the principle, the quadratic equation has a strictly negative real root iff ψ 1 > 0, ψ 2 > 0, and ψ 1 ψ 2 > 0.

3.2.6
Global stability of DFEP Theorem 3.7 The equilibrium point E 0 = (X * , 0) of system (14) is globally asymptotically stable if R l0 < 1 and conditions (i) and (ii) from the theorem in [26] are satisfied.
Proof To investigate the global stability of DFEP, we use the technique implemented by Castillo-Chavez and Song [26]. From equation (14), we can get F(X, Z) and G(X, Z), where X = (S, R l , S vl ) ∈ R 3 denotes uninfected populations and Z = (I l , I vl ) ∈ R 2 denotes infected populations: Now consider the reduced system dX dt = F(X, 0): X * = ( μ , 0, μ vl ) is a globally asymptotically stable equilibrium point for the reduced system dX dt = F(X, 0). This can be verified. From the solution of the first equation in equation (21) we obtain S(t) = μ + (S(0)μ )e -μt , which approaches μ as t − → ∞; and from the second equation of equation (21) we get S vl = μ vl + (S vl (0)μ vl )e -μ vl t , which approaches μ vd and μ vl as t − → ∞. We note that this asymptomatic dynamics is independent of the initial conditions in , therefore the convergence of the solutions of the reduced system (21) is global in . Now we compute Then G(X, Z) can be written as and we want to showĜ(X; Z) ≥ 0, which is obtained as follows: Here, μ ≥ S and μ vl ≥ S vl . Hence it is clear that equationĜ(X, Z) ≥ 0 for all (X, Z) ∈ .
Therefore, this proves that DFEP is globally asymptotically stable when 0l < 1. Hence the endemic equilibrium point does not coexist with the disease-free equilibrium point when R 0l < 1. This implies that the model exhibits forward bifurcation at R 0l = 1. From this we can say that the disease can be eliminated from the population irrespective of the initial infectious population.

Endemic equilibrium (EEP)
The endemic equilibrium is denoted by E * p = (S * , I * l , R * l , S * vl , I * vl ), and it occurs when the disease persists in the community.

Lemma 3.8 A unique endemic equilibrium exists if one of the following holds true:
( Proof To obtain it, we equate all the model equations (14) to zero. Then we obtain When we substitute S * in the first equation of the model in (14), we get where A = α 1 ρ 2 (μ + δ 2 + γ 2 )(σ 2 + μ)σ 2 γ 2 , It is important to note that the coefficient A is positive with B and C having different signs. There is unique EE whenever condition (i) or (ii) holds true. Also the other condition we did not mention is that there are two endemic equilibria if C > 0, B < 0, and B 2 -4AC > 0. However, this does not happen because we proved that DEFP of leptospirosis only model is globally asymptotically stable, which means that there is no way that the endemic equilibrium point co-exists with DFEP when R 0l < 1.

Dengue-leptospirosis co-infection model
The model in equation (1)

Invariant region
In this section, we obtain a region in which the solution of (1) is bounded. For this model the total human population is N = S + I d + I l + I dl + R d + R l + R dl . Then, after differentiating N with respect to time and substituting each equation from equation (1) to the respective place, we get If we do not consider death from leptospirosis, dengue, and the co-infection, then equation (25) becomes After solving equation (26) and evaluating it as time tends to infinity, we get Similarly, for vector population of dengue, if there is no discharge of virus from the infectious, then Also, for vector population of leptospirosis, if there is no discharge of bacteria from the infectious, then Therefore, the feasible solution set for the dengue-leptospirosis co-infection model is given by Therefore, all the solution set of equation (1) is Proof We let τ = sup{t > 0 : S 0 (t 1 ) > 0, I d0 (t 1 ) > 0, I l0 (t 1 ) > 0, I dl0 (t 1 ) > 0, R d0 (t 1 ) > 0, R l0 (t 1 ) > 0, R dl0 (t 1 ) > 0, S vd0 (t 1 ) > 0, I vd0 (t 1 ) > 0, S vl0 (t 1 ) > 0, I vl0 (t 1 ) > 0 for all t 1 ∈ [0, t]}. Since S 0 ≥ 0,

Positivity of solutions
To show that solution of the model is positive, first we take dS dt from Eq. (1).

Local stability of DFEP
Theorem 3.10 DFEP is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.
Proof The Jacobian matrix of the model at DFEP is obtained as follows: From equation (34), we can get the following characteristic polynomial: Therefore, DFEP is locally asymptotically stable iff R 0 = max{R 0d , R 0l } < 1.

Global stability of DFE Theorem 3.11
The equilibrium point E 0 = (X * , 0) of system (1) is globally asymptotically stable if R 0 ≤ 1 and conditions (i) and (ii) in the theorem [26] are satisfied.
Proof We use the technique implemented by Castillo-Chavez and Song [26]. From equation (1), we can get F(X, Z) and G(X, Z), where X = (S, R d , R l , R dl , S vd , S vl ) ∈ R 6 denotes uninfected populations and Z = (I d , I l , I dl , I vd , I vl ) ∈ R 5 denotes infected populations.
Now consider the reduced system dX dt = F(X, 0): is a globally asymptotically stable equilibrium point for the reduced system dX dt = F(X, 0). This can be verified. From the solution of the first equation in equation (35) we obtain S(t) = μ + (S(0)μ )e -μt , which approaches μ as t − → ∞; and from the second and third equation of equation (35) we get S vd = μ vd + (S vd (0)μ vd )e -μ vd t and S vl = μ vl + (S vl (0)μ vl )e -μ vl t , which approaches μ vd and μ vl as t − → ∞. We note that this asymptomatic dynamics is independent of the initial conditions in , therefore the convergence of the solutions of the reduced system (35) is global in . Now we compute Then G(X, Z) can be written as and we want to showĜ(X; Z) ≥ 0, which is obtained as follows: In equation (36),Ĝ 3 (X, Z) < 0, which leads toĜ(X, Z) < 0 for all (X, Z) ∈ . Therefore, this proves that DFEP may not be globally asymptotically stable when R 0 < 1. Hence the endemic equilibrium point coexists with DFEP when R 0 < 1. This implies that the model exhibits backward bifurcation at R 0 = 1.

Sensitivity analysis
Sensitivity is performed to identify the most dominant parameters for the spreading out as well as control of infection in the community. To go through sensitivity analysis, we follow the technique described in [27]. The sensitivity index of R 0 with respect to a parameter, say y, is given by Since R 0 = max{R 0d , R 0l }, we obtain the sensitivity analysis of R 0d and R 0l separately in the following way: The above computation of sensitivity analysis is summarized in Table 2.
From Table 2, we understand that the parameters with positive sensitivity indices, particularly β 1 , ρ 1 , c, α 1 , α 2 , and ρ 2 , have great potential in expanding dengue, leptospirosis, and their co-infection in the community. However, the parameters with negative indices, particularly δ 1 , γ 1 , δ 2 , γ 2 , μ, μ vd , and μ vl , have a great contribution in controlling the expansion of dengue, leptospirosis, and their co-infection in the community if their values are increased. From this policy makers and stakeholders are expected to act accordingly in combating dengue,leptospirosis, and their co-infection in the community.

Numerical simulations
Analytical results cannot be complete without numerical result verification. In this section, the full dengue-leptospirosis co-infection model numerical simulation is performed using Maple 18. The simulation is used for checking the effect of some parameters in the expansion as well as the control of dengue only, leptospirosis only, and co-infection of dengue and leptospirosis. For simulation purpose the parameter values in Table 3 are used.

Effect of recovery rate (γ 1 ) on dengue infectious population
In Fig. 2, we experimented with the effect of γ 1 in reducing dengue-only infectious population by maintaining the contact rate (β 1 ) constant. The figure reflects that when the values of γ 1 increase, the number of dengue only infectious population is going down. From this we should concentrate on improving recovery rates either by treating infected populations or by raising individuals' immunity levels to dengue disease in the population. It should be viewed by policy makers as a mitigation strategy.

Effect of recovery rate (γ 2 ) on leptospirosis infectious population
In Fig. 3, it is shown that γ 2 plays a significant role in bringing down the leptospirosis infection. When the value of γ 2 increased from 0.017 to 0.087, the amount of infectious population due to leptospirosis decreased, where the contact rate is kept constant, which is α 1 = α 2 . It can also be used by policy makers as a tool for mitigation.

Effect of dengue contact rate (β 1 ) on co-infectious population
As it is shown in Fig. 4, the contact rate of leptospirosis (α 1 = α 2 ) and the recovery rate of co-infectious population (ε) are kept constant. The figure reflects that as the value of contact rate of dengue increased, the co-infectious population increased, which means increased expansion of co-infection. From this we can see that decreasing of the dengue contact rate is significant in the controlling of co-infection transmission. Therefore, stakeholders are expected to work on decreasing the contact rate of susceptible humans and dengue vector by using either bednet or chemical or using an appropriate method of prevention mechanism to bring down the expansion of co-infection in the community.

Effect of leptospirosis contact rate (α 1 and α 2 ) on co-infectious population
Similarly, we investigated the effect of leptospirosis contact rate (α 1 = α 2 ) in the expansion of dengue-leptospirosis co-infection while keeping the recovery rate of co-infection (ε) constant. Figure 5 shows that co-infectious population decreases as the leptospirosis contact rate is decreasing, by keeping dengue contact rate (β 1 ) and (ε) constant. This implies that, in order to mitigate the co-infection, it is advisable to reduce the rate of contact with infectious humans, and that the host vector of leptospirosis is crucial.

Effect of recovery rate of dengue-leptospirosis (ε) on co-infectious population
Here, we experimented on the effect of recovery rate of dengue and leptospirosis (ε) on the co-infectious population. As we explained in the model description, due to treatment or other mechanisms, co-infectious population either recover from dengue only or from leptospirosis only or from both diseases with their own probability and join their respective recovered compartment. Therefore, Fig. 6 shows that increasing the rate of recovery of the co-infectious population has a great advantage in reducing both diseases in the population.

Figure 5
Effect of leptospirosis contact rate on co-infectious population Figure 6 Effect of recovery rate of dengue-leptospirosis on co-infectious population

Conclusions
The deterministic co-infection model for dengue-leptospirosis disease was developed using ordinary differential equations, and the population is subdivided into eleven compartments. The qualitative analysis of the model was done by splitting the full model into two, which are dengue only and leptospirosis only models. The analysis of the model shows that there exists a region where the model is mathematically and epidemiologically well posed. Basic reproduction numbers, disease-free equilibrium, endemic equilibria, stability analysis of equilibrium points, and sensitivity analysis of basic reproduction of dengue only, leptospirosis only, and the full model were analyzed in their respective order. Numerically, we experimented on the effect of basic parameters in the expansion or control of dengue only, leptospirosis only, and co-infectious diseases. From the result, we conclude that an increase in the rate of dengue recovery contributes greatly to reducing dengue infection in the community. Similarly, increasing the recovery rate for leptospirosis also contributes to the reduction of leptospirosis. The rate of recovery for co-infection also has an effect on reducing co-infectious population if its value has been increased. The other finding obtained in this section is that decreasing either dengue or leptospirosis contact rate has a great influence on controlling co-infection of dengue and leptospirosis in the population.