# Transmission Dynamics of Lymphatic Filariasis: A Mathematical Approach

Authors: C. P. Bhunu; S. Mushayabasa
Journal: ISRN Biomathematics
 (2012)

---

## Abstract

An epidemiological model for the spread of lymphatic filariasis, a mosquito-borne infection, is developed and analysed. The epidemic thresholds known as the reproduction number and equilibria for the model are determined and stabilities analysed. Results from the analysis of the reproduction number suggest that treatment will somehow contribute to a reduction in lymphatic filariasis cases, but what it does not show is the magnitude of the reduction, a part answered by the numerical simulations. Numerical simulations show that even when all lymphatic filariasis cases displaying elephantiasis
symptoms are put on treatment it will not be able to eradicate the disease. This result suggests that effective control of lymphatic filariasis may lie in treatment for those displaying symptoms as well as chemoprophylaxis for the exposed.


---

## Body

## 1. Introduction

Lymphatic filariasis, a debilitating disease, is one of the most prevalent and yet one of the most neglected tropical diseases with serious economic and social consequences [1, 2]. Lymphatic filariasis affects women, men and children of all ages. It is a mosquito-borne disease caused by tissue-dwelling nematodes of Brugia malayi, Brugia timori, and Wuchereria bancrofti species [1, 3] and is estimated to affect about 120 million people worldwide [4–6].Wuchereria bancrofti is responsible for 90% of the cases and is found throughout the tropical and subtropical areas of the world; Brugia malayi is confined to southeast and eastern Asia; Brugia timori is found only in Timor and its adjacent islands [7]. Infection leads to lymphedema, a buildup of fluid due to impaired function of the lymph vessels, in only a small proportion people, even in areas of intense transmission [8], as most people with long-term infections are clinically asymptomatic. Recurrent bacterial infections in some lymphedema patients lead to elephantiasis [9].Filarial parasites are a major cause of morbidity and therefore hinder socioeconomic growth in parts of Asia, Africa, and the Western Pacific [1, 10]. Despite improved knowledge of pathology of lymphatic filariasis and existence of the drugs diethylcarbamazine and albendazole necessary to treat lymphatic filariasis, it continues to be a major public health problem in tropical and subtropical countries. Lymphatic filariasis is more common in regions that have a higher incidence of poverty [11] making it a disease of the poor and serves as an indicator of underdevelopment [1]. Surveys in Ghana have indicated that bancroftian filariasis is present in most parts of the country with considerable variations in prevalence [12]. In countries where lymphatic filariasis is well established, the prevalence of infection continues to increase due to unplanned growth of cities and water resource development such as irrigation, which creates numerous breeding sites for the mosquitoes that transmit the disease [12].A number of mathematical models have looked into malaria a mosquito-borne infection [13–19], to mention just a few, but only a few have looked into lymphatic filariasis [20–24]. Here the authors propose and analyse the transmission dynamics of lymphatic filariasis using differential equations to explore if treatment for those symptoms alone will be able to keep the infection under control.


## 2. Model Description

The model subdivides the human population based on lymphatic filariasis status into the following subpopulations: susceptible humansSh(t), infected and not showing signs of elephantiasis (exposed) Eh(t), infected and displaying elephantiasis symptoms Ih(t). Thus, the total human population is given by
(1)Nh(t)=Sh(t)+Eh(t)+Ih(t).
The mosquito population is divided into the following subgroups: susceptible mosquitoes Sm(t) and infected mosquitoes Im(t), so the total mosquito population is given by
(2)Nm(t)=Sm(t)+Im(t).
The mosquitoes and human beings are recruited into their susceptible corresponding populations at rates Λm and Λh, respectively. Mosquitoes experience natural death rate at a rate μm which is proportional to the number in each mosquito class. Similarly, human beings experience natural death at a rate μh, which is proportional to the number in each human class.The mosquito ingests microfilariae when biting a human who is infected with filariasis (elephantiasis causing nematodes) at a rate(3)λh(t)=βh[θhEh(t)+Ih(t)]Nh(t).
Here, βh is the average number of mosquito bites which cause transmission of disease from infected human to susceptible mosquito; θh∈(0,1) accounts for reduced number of microfilariae in the blood stream of individuals infected but not showing elephantiasis symptoms. Upon getting infected, susceptible mosquitoes enter the infected class Im(t). Microfilariae pass through mosquito gut into hemocoel and develop into filariform juveniles. Filariform juveniles escape from mosquito’s proboscis when the insect is feeding and then penetrate wound structure of a human being at a rate λm(t) where
(4)λm(t)=βmIm(t)Nm(t),
where βm is the average number of mosquito bites which cause transmission of disease from infectious mosquito to susceptible human per mosquito. Thus, humans are infected at a rate λm(t) following a bite by a mosquito to move into the exposed class Eh(t) (infected and not showing signs of filariasis). Individuals in Eh(t) progress to the stage of showing filarasis symptoms Ih(t) at rate ρ. However, some progress to the Ih(t) as a result of reinfection at a rate λm(t). Individuals in Ih(t)-class are treated using Diethycarbamzime (DEC) and albendazole drugs at a rate α to move back into the susceptible class Sh(t), since recovery from filarasis is not permanent. With filarasis, there is no disease-induced death. The structure of model is given in Figure 1.Figure 1
Structure of the model.Unless stated otherwise, values for the parameters in the simulations are given in Table1 which are values used in malaria models which is also spread by the same vector as data on this neglected tropical infection is hard to come by.Table 1
Model parameters and their interpretations.
Definition
Symbol
Units
Point estimate
Range
Source
Mosquito recruitment rate
Λ
m
Mosquitoes/day
4.227
·
1
0
4
—
a*
Mosquito mortality rate
μ
m
1/day
3.623
0.01–5.0
a*
Average number of mosquito bites which cause transmission of disease from infected human to susceptible mosquito
β
h
1/day
3.576
1.52–5
a*
Modification parameter
θ
h
—
0.25
0.125–1
Assumed
Human recruitment rate
Λ
h
Humans/day
0.0024
·
2
·
1
0
4
—
a*
Human mortality rate
μ
h
1/day
0.00048
0–0.25
a*
Average number of mosquito bites which cause transmission of disease from infectious mosquito to susceptible human per mosquito
β
m
1/day
0.091
0.01–0.75
a*
Treatment rate
α
1/day
0.125
0.1–0.75
Assumed
Per capita rate of progression of human from the exposed class to the infectious class
ρ
1/day
0.01
0–0.25
Assumed
a* denotes parameter values from Labadin et al., [25].Due to lack of data to calibrate the model and for parameter estimation, other parameter values are assumed within realistic ranges for illustrative purpose. Note that that “Range” refers to values used in the sensitivity analysis. Based on the given assumptions, the following system of differential equations describes the model:(5)Sm′(t)=Λm-(λh+μm)Sm,Im′(t)=λhSm-μmIm,Sh′(t)=Λh-(λm+μh)Sh+αIh,Eh′(t)=λmSh-(λm+ρ+μh)Eh,Ih′(t)=(λm+ρ)Eh-(μh+α)Ih.
All parameters and state variables for model system  are assumed to be nonnegative to be consistent with human and mosquito populations. All feasible solutions of model system  are positive and eventually enter the invariant attracting region
(6)Ω={(Sm,Im)∈ℝ+2:Nm≤Λmμm,(Sh,Eh,Ih)∈ℝ+3:Nh≤Λhμh.
It is sufficient to consider solutions in Ω. Existence, uniqueness, and continuation results for system  hold in this region, and all solutions starting in Ω remain there for all t≥0. Hence,  is mathematically and epidemiologically well posed, and it is sufficient to consider the dynamics of the flow generated by the model system  in Ω.


### 2.1. Disease-Free Equilibrium and Stability Analysis

The disease-free equilibrium of model system is given by
(7)ℰ0=(Sm0,Im0,Sh0,Eh0,Ih0)=(Λmμm,0,Λhμh,0,0).
Following van den Driessche and Watmough [26], we have
(8)ℛE=βmβh(θh(α+μh)+ρ)μm(α+μh)(ρ+μh),
which is the number secondary lymphatic filarasis cases caused by one infected individual/mosquito in a completely susceptible population. Local stability of the disease-free equilibrium is assured by Theorem 2 [26].Theorem 1.
The disease-free equilibriumℰ0 for model system  is locally asymptotically stable if ℛE<1 and unstable otherwise.


#### 2.1.1. Analysis of the Reproduction NumberℛE

ℛ
E is a decreasing function of μm as limμm→∞ℛE=0, suggesting that mechanisms which increase mortality of mosquitoes like spraying will have a positive impact on controlling lymphatic filariasis. Now let us explore the impact of treating elephantiasis cases only by partially differentiating ℛE with respect to α to obtain
(9)∂ℛE∂α=-βmβhρ2ℛEμm(μh+ρ)(μh+α)2<0.
The fact that ∂ℛE/∂α<0 suggests that increase in treatment of elephantiasis cases results in a decrease of the lymphatic filariasis infections in the population.


### 2.2. Existence and Stability Analysis of the Endemic Equilibrium

The endemic equilibrium is given byℰ*=(Sm*,Im*,Sh*,Eh*,Ih*), where
(10)Sm*=Λmμm+λh*,Im*=Λmλh*μm(μm+λh*),Nm*=Λmμm,Sh*=Λh(λm*+μh+ρ)(μh+α)μh[λm*2+λm*(2(μh+α)+ρ)+(μh+ρ)(μh+α)],Eh*=Λhλm*(μh+α)μh[λm*2+λm*(2(μh+α)+ρ)+(μh+ρ)(μh+α)],Ih*=Λhλm*(λm*+ρ)μh[λm*2+λm*(2(μh+α)+ρ)+(μh+ρ)(μh+α)],Nh*=Λhμh.
Substituting  into the force of infection λm* in (4), we obtain
(11)λm*=βmIm*Nm*=βmλh*μm+λh*⇒λh*=μmλm*βm-λm*.
Also, substituting  into the force of infection λh* in (3), we obtain
(12)λh*=βh[λm*2+λm*(θh(μh+α)+ρ)]λm*2+λm*(2(μh+α)+ρ)+(μh+α)(μh+ρ).
Substituting (11) into (12), we have λm*h(λm*)=0, with λm*=0 corresponding to the disease-free equilibrium and h(λm*)=0 corresponding to the endemic equilibrium with
(13)h(λm*)=Aλm*2+Bλm*+C=0,A=μm+βhμm(μh+α)(μh+ρ),C=1-ℛE2,B=(μh+α)(2μm+βhθh)+ρ(μm+βh)-βhβmμm(μh+α)(μh+ρ).
Clearly A is positive throughout, C>0(C<0) only when ℛE<1(ℛE>1) and B can be positive or negative depending on whether (μh+α)(2μm+βhθh)+ρ(μm+βh)> or <βhβm. Solving h(λm*)=0, we obtain
(14)λm*=-B±B2-4AC2A.
The analysis of (14) is Theorem 2.Theorem 2.
Model system has  (i)
a unique endemic equilibrium ifC<0(ℛE>1);(ii)
a backward bifurcation ifB<0, C>0(ℛE<1).We now employ the Centre Manifold Theory [27] to analyse the stability of this equilibrium point as described in Theorem 4.1 [28], to establish the local asymptotic stability of the endemic equilibrium. Let us make the following change of variables in order to apply the Center Manifold Theory Sm=x1,Im=x2, so that Nm=∑n=12xn and Sh=x3,Eh=x4,Ih=x5, so that Nh=∑n=13x2+n. We now use the vector notation X=(x1,x2,x3,x4,x5)T. Then, model system  can be written in the form dX/dt=F=(f1,f2,f3,f4,f5)T, such that
(15)x1′(t)=f1=Λm-βh(θhx4+x5)∑n=13xn+2x1-μmx1,x2′(t)=f2=βh(θhx4+x5)∑n=13xn+2x1-μmx2,x3′(t)=f3=Λh-βmx2∑n=12xnx3-μhx3+αx5,x4′(t)=f4=βmx2∑n=12xnx3-βmx2∑n=12xnx4-(ρ+μh)x4,x5′(t)=f5=βmx2∑n=12xnx4+ρx4-(μh+α)x5.
The Jacobian matrix of system  at ℰ0 is given by
(16)J(ℰ0)=[-μm00-βhθhΛmμhμmΛh-βhΛmμhμmΛh0-μm0βhθhΛmμhμmΛhβhΛmμhμmΛh0-βmΛhμmμhΛm-μh0α0βmΛhμmμhΛm0-(μh+ρ)0000ρ-(μh+α)],
from which it can be shown that the reproduction number is
(17)ℛE=βmβh(θh(α+μh)+ρ)μm(α+μh)(ρ+μh).
Now let us consider βm=κβh, regardless of whether κ∈(0,1) or κ≥1. Taking βh as a bifurcation parameter and if we consider the case ℛE=1 and solve for βh, we obtain
(18)βh=β*=μm(μh+α)(μh+ρ)κ(θh(μh+α)+ρ).
Note that the linearized system of the transformed equation  with the bifurcation point β* has a simple zero eigenvalue. Hence, the Centre Manifold Theory [27] can be used to analyze the dynamics of  near βh=β*. It can be shown that the Jacobian of  at βh=β* has a right eigenvector associated with the zero eigenvalue given by u=[u1,u2,u3,u4,u5]T, where
(19)u1=-β*Λmμh(θh(μh+α)+ρ)μm2Λh(μh+α)u4,u2=μhΛm(μh+ρ)Λhμmκβ*u4,u3=-(μh+ρ+α)u4,u4>0,u5=ρμh+αu4.
The left eigenvector of J(ℰ0) associated with the zero eigenvalue at βh=β* is given by v=[v1,v2,v3,v4,v5]T, where
(20)v1=0,v2=κβ*ΛhμhΛmv4>0,v3=0,v4>0,v5=κβ*2μm(μh+α)v4>0.


#### 2.2.1. Computation of Bifurcation Parametersa and b

The sign ofb is associated with the following nonvanishing partial derivatives of F(21)∂2f2∂β*∂x4=ΛmμhθhμmΛh,∂2f2∂β*∂x5=ΛmμhμmΛh,∂2f4∂β*∂x2=κΛhμmμhΛm.
Thus, it follows from , (20), and (21) that
(22)b=∑i,j=15‍viuj∂2fi∂β*∂xj=(κβ*(θh(μh+α)+ρ)μm(μh+α)+μh+ρβ*)u4v4>0.
For system , the associated nonzero partial derivatives of F: at the disease-free equilibrium associated with a are given by
(23)∂2f2∂x1∂x4=∂2f2∂x4∂x1=β*θhμhΛh,∂2f2∂x1∂x5=∂2f2∂x5∂x1=β*μhΛh,∂2f2∂x3∂x4=∂2f2∂x4∂x3=-β*θhμh2ΛmΛh2μm,∂2f2∂x3∂x5=∂2f2∂x5∂x3=-β*μh2ΛmΛh2μm,∂2f2∂x42=-2β*θhμh2ΛmΛh2μm,∂2f2∂x4∂x5=∂2f2∂x5∂x4=-β*(θh+1)μh2ΛmΛh2μm,∂2f2∂x52=-2β*μh2ΛmΛh2μm,∂2f4∂x1∂x2=∂2f4∂x2∂x1=-κβ*μm2ΛhΛm2μh,∂2f4∂x22=-2κβ*μm2ΛhΛm2μh,∂2f4∂x2∂x3=∂2f4∂x3∂x2=κβ*μmΛm,∂2f4∂x2∂x4=∂2f4∂x4∂x2=-κβ*μmΛm,∂2f5∂x2∂x4=∂2f5∂x4∂x2=κβ*μmΛm.
It follows from , (20), and  that(24)a=∑i,j,k=15‍viujuk∂2fi∂xj∂xk=κβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2]-[μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)]κβ*Λhγ12μm2u42v4,with γ1=μh+α,γ2=μh+ρ,γ3=μh+α+ρ,γ4=θh(μh+α)+ρ. It follows that a> or <0 if κβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2]> or <[μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)].Thus,b>0 and a>0 or a<0 depending on whether κβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2] is greater or smaller than [μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)]. This sign of b may be expected in general for epidemic models because, in essence, using β* as a bifurcation parameter often ensures b>0. Using Theorem 4.1 [28] items (i) and (iv), we establish Theorem 3.Theorem 3.
Ifκβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2]>[μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)], a>0, then model system  undergoes a backward bifurcation at ℛE=1, otherwise a<0 and a unique endemic equilibrium ℰ* guaranteed by Theorem 4.1 [28] is locally asymptotically stable for ℛE>1, but close to 1.The phenomenon of backward bifurcation in disease models, where a stable endemic equilibrium coexists with a stable disease-free equilibrium when the associated reproduction number is less than unity, has important implications for disease control [29]. In such a case, having the reproduction number being less than unity becomes only a necessary but are not sufficient condition for disease elimination. Backward bifurcation in this case is caused by reinfection of the exposed individuals. In the absence of reinfection model system (7.3) will have a unique endemic equilibrium which exists whenever ℛE>1(a<0). Thus, for a<0, the model exhibits a forward bifurcation. The bifurcations which occur for different signs of a are shown in Figure 2 using parameters values in Table 1.Figure 2
Bifurcation diagram.


## 2.1. Disease-Free Equilibrium and Stability Analysis

The disease-free equilibrium of model system is given by
(7)ℰ0=(Sm0,Im0,Sh0,Eh0,Ih0)=(Λmμm,0,Λhμh,0,0).
Following van den Driessche and Watmough [26], we have
(8)ℛE=βmβh(θh(α+μh)+ρ)μm(α+μh)(ρ+μh),
which is the number secondary lymphatic filarasis cases caused by one infected individual/mosquito in a completely susceptible population. Local stability of the disease-free equilibrium is assured by Theorem 2 [26].Theorem 1.
The disease-free equilibriumℰ0 for model system  is locally asymptotically stable if ℛE<1 and unstable otherwise.


### 2.1.1. Analysis of the Reproduction NumberℛE

ℛ
E is a decreasing function of μm as limμm→∞ℛE=0, suggesting that mechanisms which increase mortality of mosquitoes like spraying will have a positive impact on controlling lymphatic filariasis. Now let us explore the impact of treating elephantiasis cases only by partially differentiating ℛE with respect to α to obtain
(9)∂ℛE∂α=-βmβhρ2ℛEμm(μh+ρ)(μh+α)2<0.
The fact that ∂ℛE/∂α<0 suggests that increase in treatment of elephantiasis cases results in a decrease of the lymphatic filariasis infections in the population.


## 2.1.1. Analysis of the Reproduction NumberℛE

ℛ
E is a decreasing function of μm as limμm→∞ℛE=0, suggesting that mechanisms which increase mortality of mosquitoes like spraying will have a positive impact on controlling lymphatic filariasis. Now let us explore the impact of treating elephantiasis cases only by partially differentiating ℛE with respect to α to obtain
(9)∂ℛE∂α=-βmβhρ2ℛEμm(μh+ρ)(μh+α)2<0.
The fact that ∂ℛE/∂α<0 suggests that increase in treatment of elephantiasis cases results in a decrease of the lymphatic filariasis infections in the population.


## 2.2. Existence and Stability Analysis of the Endemic Equilibrium

The endemic equilibrium is given byℰ*=(Sm*,Im*,Sh*,Eh*,Ih*), where
(10)Sm*=Λmμm+λh*,Im*=Λmλh*μm(μm+λh*),Nm*=Λmμm,Sh*=Λh(λm*+μh+ρ)(μh+α)μh[λm*2+λm*(2(μh+α)+ρ)+(μh+ρ)(μh+α)],Eh*=Λhλm*(μh+α)μh[λm*2+λm*(2(μh+α)+ρ)+(μh+ρ)(μh+α)],Ih*=Λhλm*(λm*+ρ)μh[λm*2+λm*(2(μh+α)+ρ)+(μh+ρ)(μh+α)],Nh*=Λhμh.
Substituting  into the force of infection λm* in (4), we obtain
(11)λm*=βmIm*Nm*=βmλh*μm+λh*⇒λh*=μmλm*βm-λm*.
Also, substituting  into the force of infection λh* in (3), we obtain
(12)λh*=βh[λm*2+λm*(θh(μh+α)+ρ)]λm*2+λm*(2(μh+α)+ρ)+(μh+α)(μh+ρ).
Substituting (11) into (12), we have λm*h(λm*)=0, with λm*=0 corresponding to the disease-free equilibrium and h(λm*)=0 corresponding to the endemic equilibrium with
(13)h(λm*)=Aλm*2+Bλm*+C=0,A=μm+βhμm(μh+α)(μh+ρ),C=1-ℛE2,B=(μh+α)(2μm+βhθh)+ρ(μm+βh)-βhβmμm(μh+α)(μh+ρ).
Clearly A is positive throughout, C>0(C<0) only when ℛE<1(ℛE>1) and B can be positive or negative depending on whether (μh+α)(2μm+βhθh)+ρ(μm+βh)> or <βhβm. Solving h(λm*)=0, we obtain
(14)λm*=-B±B2-4AC2A.
The analysis of (14) is Theorem 2.Theorem 2.
Model system has  (i)
a unique endemic equilibrium ifC<0(ℛE>1);(ii)
a backward bifurcation ifB<0, C>0(ℛE<1).We now employ the Centre Manifold Theory [27] to analyse the stability of this equilibrium point as described in Theorem 4.1 [28], to establish the local asymptotic stability of the endemic equilibrium. Let us make the following change of variables in order to apply the Center Manifold Theory Sm=x1,Im=x2, so that Nm=∑n=12xn and Sh=x3,Eh=x4,Ih=x5, so that Nh=∑n=13x2+n. We now use the vector notation X=(x1,x2,x3,x4,x5)T. Then, model system  can be written in the form dX/dt=F=(f1,f2,f3,f4,f5)T, such that
(15)x1′(t)=f1=Λm-βh(θhx4+x5)∑n=13xn+2x1-μmx1,x2′(t)=f2=βh(θhx4+x5)∑n=13xn+2x1-μmx2,x3′(t)=f3=Λh-βmx2∑n=12xnx3-μhx3+αx5,x4′(t)=f4=βmx2∑n=12xnx3-βmx2∑n=12xnx4-(ρ+μh)x4,x5′(t)=f5=βmx2∑n=12xnx4+ρx4-(μh+α)x5.
The Jacobian matrix of system  at ℰ0 is given by
(16)J(ℰ0)=[-μm00-βhθhΛmμhμmΛh-βhΛmμhμmΛh0-μm0βhθhΛmμhμmΛhβhΛmμhμmΛh0-βmΛhμmμhΛm-μh0α0βmΛhμmμhΛm0-(μh+ρ)0000ρ-(μh+α)],
from which it can be shown that the reproduction number is
(17)ℛE=βmβh(θh(α+μh)+ρ)μm(α+μh)(ρ+μh).
Now let us consider βm=κβh, regardless of whether κ∈(0,1) or κ≥1. Taking βh as a bifurcation parameter and if we consider the case ℛE=1 and solve for βh, we obtain
(18)βh=β*=μm(μh+α)(μh+ρ)κ(θh(μh+α)+ρ).
Note that the linearized system of the transformed equation  with the bifurcation point β* has a simple zero eigenvalue. Hence, the Centre Manifold Theory [27] can be used to analyze the dynamics of  near βh=β*. It can be shown that the Jacobian of  at βh=β* has a right eigenvector associated with the zero eigenvalue given by u=[u1,u2,u3,u4,u5]T, where
(19)u1=-β*Λmμh(θh(μh+α)+ρ)μm2Λh(μh+α)u4,u2=μhΛm(μh+ρ)Λhμmκβ*u4,u3=-(μh+ρ+α)u4,u4>0,u5=ρμh+αu4.
The left eigenvector of J(ℰ0) associated with the zero eigenvalue at βh=β* is given by v=[v1,v2,v3,v4,v5]T, where
(20)v1=0,v2=κβ*ΛhμhΛmv4>0,v3=0,v4>0,v5=κβ*2μm(μh+α)v4>0.


### 2.2.1. Computation of Bifurcation Parametersa and b

The sign ofb is associated with the following nonvanishing partial derivatives of F(21)∂2f2∂β*∂x4=ΛmμhθhμmΛh,∂2f2∂β*∂x5=ΛmμhμmΛh,∂2f4∂β*∂x2=κΛhμmμhΛm.
Thus, it follows from , (20), and (21) that
(22)b=∑i,j=15‍viuj∂2fi∂β*∂xj=(κβ*(θh(μh+α)+ρ)μm(μh+α)+μh+ρβ*)u4v4>0.
For system , the associated nonzero partial derivatives of F: at the disease-free equilibrium associated with a are given by
(23)∂2f2∂x1∂x4=∂2f2∂x4∂x1=β*θhμhΛh,∂2f2∂x1∂x5=∂2f2∂x5∂x1=β*μhΛh,∂2f2∂x3∂x4=∂2f2∂x4∂x3=-β*θhμh2ΛmΛh2μm,∂2f2∂x3∂x5=∂2f2∂x5∂x3=-β*μh2ΛmΛh2μm,∂2f2∂x42=-2β*θhμh2ΛmΛh2μm,∂2f2∂x4∂x5=∂2f2∂x5∂x4=-β*(θh+1)μh2ΛmΛh2μm,∂2f2∂x52=-2β*μh2ΛmΛh2μm,∂2f4∂x1∂x2=∂2f4∂x2∂x1=-κβ*μm2ΛhΛm2μh,∂2f4∂x22=-2κβ*μm2ΛhΛm2μh,∂2f4∂x2∂x3=∂2f4∂x3∂x2=κβ*μmΛm,∂2f4∂x2∂x4=∂2f4∂x4∂x2=-κβ*μmΛm,∂2f5∂x2∂x4=∂2f5∂x4∂x2=κβ*μmΛm.
It follows from , (20), and  that(24)a=∑i,j,k=15‍viujuk∂2fi∂xj∂xk=κβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2]-[μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)]κβ*Λhγ12μm2u42v4,with γ1=μh+α,γ2=μh+ρ,γ3=μh+α+ρ,γ4=θh(μh+α)+ρ. It follows that a> or <0 if κβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2]> or <[μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)].Thus,b>0 and a>0 or a<0 depending on whether κβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2] is greater or smaller than [μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)]. This sign of b may be expected in general for epidemic models because, in essence, using β* as a bifurcation parameter often ensures b>0. Using Theorem 4.1 [28] items (i) and (iv), we establish Theorem 3.Theorem 3.
Ifκβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2]>[μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)], a>0, then model system  undergoes a backward bifurcation at ℛE=1, otherwise a<0 and a unique endemic equilibrium ℰ* guaranteed by Theorem 4.1 [28] is locally asymptotically stable for ℛE>1, but close to 1.The phenomenon of backward bifurcation in disease models, where a stable endemic equilibrium coexists with a stable disease-free equilibrium when the associated reproduction number is less than unity, has important implications for disease control [29]. In such a case, having the reproduction number being less than unity becomes only a necessary but are not sufficient condition for disease elimination. Backward bifurcation in this case is caused by reinfection of the exposed individuals. In the absence of reinfection model system (7.3) will have a unique endemic equilibrium which exists whenever ℛE>1(a<0). Thus, for a<0, the model exhibits a forward bifurcation. The bifurcations which occur for different signs of a are shown in Figure 2 using parameters values in Table 1.Figure 2
Bifurcation diagram.


## 2.2.1. Computation of Bifurcation Parametersa and b

The sign ofb is associated with the following nonvanishing partial derivatives of F(21)∂2f2∂β*∂x4=ΛmμhθhμmΛh,∂2f2∂β*∂x5=ΛmμhμmΛh,∂2f4∂β*∂x2=κΛhμmμhΛm.
Thus, it follows from , (20), and (21) that
(22)b=∑i,j=15‍viuj∂2fi∂β*∂xj=(κβ*(θh(μh+α)+ρ)μm(μh+α)+μh+ρβ*)u4v4>0.
For system , the associated nonzero partial derivatives of F: at the disease-free equilibrium associated with a are given by
(23)∂2f2∂x1∂x4=∂2f2∂x4∂x1=β*θhμhΛh,∂2f2∂x1∂x5=∂2f2∂x5∂x1=β*μhΛh,∂2f2∂x3∂x4=∂2f2∂x4∂x3=-β*θhμh2ΛmΛh2μm,∂2f2∂x3∂x5=∂2f2∂x5∂x3=-β*μh2ΛmΛh2μm,∂2f2∂x42=-2β*θhμh2ΛmΛh2μm,∂2f2∂x4∂x5=∂2f2∂x5∂x4=-β*(θh+1)μh2ΛmΛh2μm,∂2f2∂x52=-2β*μh2ΛmΛh2μm,∂2f4∂x1∂x2=∂2f4∂x2∂x1=-κβ*μm2ΛhΛm2μh,∂2f4∂x22=-2κβ*μm2ΛhΛm2μh,∂2f4∂x2∂x3=∂2f4∂x3∂x2=κβ*μmΛm,∂2f4∂x2∂x4=∂2f4∂x4∂x2=-κβ*μmΛm,∂2f5∂x2∂x4=∂2f5∂x4∂x2=κβ*μmΛm.
It follows from , (20), and  that(24)a=∑i,j,k=15‍viujuk∂2fi∂xj∂xk=κβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2]-[μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)]κβ*Λhγ12μm2u42v4,with γ1=μh+α,γ2=μh+ρ,γ3=μh+α+ρ,γ4=θh(μh+α)+ρ. It follows that a> or <0 if κβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2]> or <[μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)].Thus,b>0 and a>0 or a<0 depending on whether κβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2] is greater or smaller than [μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)]. This sign of b may be expected in general for epidemic models because, in essence, using β* as a bifurcation parameter often ensures b>0. Using Theorem 4.1 [28] items (i) and (iv), we establish Theorem 3.Theorem 3.
Ifκβ*2μmγ1[γ4(μhγ2+κβ*γ3)+μmμhκβ*γ2]>[μm2γ12μhγ2(κβ*(1+γ3)+γ2)+κ2β*2γ4(μmγ3+β*2μhγ4)], a>0, then model system  undergoes a backward bifurcation at ℛE=1, otherwise a<0 and a unique endemic equilibrium ℰ* guaranteed by Theorem 4.1 [28] is locally asymptotically stable for ℛE>1, but close to 1.The phenomenon of backward bifurcation in disease models, where a stable endemic equilibrium coexists with a stable disease-free equilibrium when the associated reproduction number is less than unity, has important implications for disease control [29]. In such a case, having the reproduction number being less than unity becomes only a necessary but are not sufficient condition for disease elimination. Backward bifurcation in this case is caused by reinfection of the exposed individuals. In the absence of reinfection model system (7.3) will have a unique endemic equilibrium which exists whenever ℛE>1(a<0). Thus, for a<0, the model exhibits a forward bifurcation. The bifurcations which occur for different signs of a are shown in Figure 2 using parameters values in Table 1.Figure 2
Bifurcation diagram.


## 3. Numerical Simulations

In this section, we make use of Matlab to analyse the effect of varying initial infected mosquito population, for different reproduction numbers of model system using model parameters in Table 1.Figure3 shows that whenever the reproduction number is greater than unity, then all the populations (mosquito and human) converge to their corresponding respective endemic equilibrium point regardless of the size of initial mosquito population size.Simulations of model system for varying initial infected mosquito population when ℛE>1(ℛE=1.809). Parameter values used are in Table 1.
(a)
(b)In Figure4, the effect of varying infected mosquito population size when all human individuals displaying signs and symptoms of elephantiasis are put on treatment is illustrated. It shows that even when the reproduction number is less than unity the individuals displaying elephantiasis symptoms will not disappear in the population. This result suggests that effective control of elephantiasis may lie in treatment of both humans in the exposed and infectious classes. Reinfection of the people in the exposed class makes it difficult to eliminate lymphatic filariasis through treatment of the infectives alone.Simulations of model system for varying initial infected mosquito population when ℛE<1(ℛE=0.715). Parameter values used are in Table 1.
(a)
(b)
(c)


### 3.1. Sensitivity Analysis

In order to investigate the effects of variations inℛE to its constituent parameters, we used Latin Hypercube Sampling and Partial Rank Correlation Coefficients (PRCCs) with 1000 simulations per run. Latin Hypercube Sampling is a statistical sampling method that allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter [30]. PRCCs illustrate the degree of the effect that each parameter has on the outcome.Figure5 illustrates the PRCCs using ℛE as an output variable. The parameter with the greatest effect on the outcome is the death rate of mosquitoes suggesting that mechanisms which increase death rate of mosquitoes will have a positive impact on lymphatic filariasis control. Surprisingly, the rate of progressing from the exposed class to the infectious class decreases ℛE more than treatment. This result turns to suggest that treatment of people displaying symptoms of elephantiasis symptoms alone is not enough to keep the disease under control.Figure 5
Partial rank correlation coefficients showing the effect of parameter variations onℛE using ranges in Table 1. Parameters with positive PRCCs will increase ℛE when they are increased, whereas parameters with negative PRCCs will decrease ℛE when they are increased.Figure6 illustrates the effect that varying six sample parameters will have on ℛE. If the mosquito death rates are sufficiently high, then the disease prevalence can be reduced (ℛE→1), but not sufficiently enough to control it (ℛE<1). Similarly, high treatment levels result in ℛE→1 but are not sufficient enough to have ℛE<1 necessary for disease control. Increases in the average number of mosquito bites (βm,βh) increase ℛE suggesting an increase in lymphatic filariasis cases. This further suggests that other mechanisms which reduce mosquito bites like spraying, use of mosquito nets, and wearing clothes which cover much of the body may be necessary to control lymphatic filariasis.Monte Carlo simulations of 1000 sample values for six illustrative parameters (βh,βm,ρ,α,andθh) chosen via Latin Hypercube Sampling.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.1. Sensitivity Analysis

In order to investigate the effects of variations inℛE to its constituent parameters, we used Latin Hypercube Sampling and Partial Rank Correlation Coefficients (PRCCs) with 1000 simulations per run. Latin Hypercube Sampling is a statistical sampling method that allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter [30]. PRCCs illustrate the degree of the effect that each parameter has on the outcome.Figure5 illustrates the PRCCs using ℛE as an output variable. The parameter with the greatest effect on the outcome is the death rate of mosquitoes suggesting that mechanisms which increase death rate of mosquitoes will have a positive impact on lymphatic filariasis control. Surprisingly, the rate of progressing from the exposed class to the infectious class decreases ℛE more than treatment. This result turns to suggest that treatment of people displaying symptoms of elephantiasis symptoms alone is not enough to keep the disease under control.Figure 5
Partial rank correlation coefficients showing the effect of parameter variations onℛE using ranges in Table 1. Parameters with positive PRCCs will increase ℛE when they are increased, whereas parameters with negative PRCCs will decrease ℛE when they are increased.Figure6 illustrates the effect that varying six sample parameters will have on ℛE. If the mosquito death rates are sufficiently high, then the disease prevalence can be reduced (ℛE→1), but not sufficiently enough to control it (ℛE<1). Similarly, high treatment levels result in ℛE→1 but are not sufficient enough to have ℛE<1 necessary for disease control. Increases in the average number of mosquito bites (βm,βh) increase ℛE suggesting an increase in lymphatic filariasis cases. This further suggests that other mechanisms which reduce mosquito bites like spraying, use of mosquito nets, and wearing clothes which cover much of the body may be necessary to control lymphatic filariasis.Monte Carlo simulations of 1000 sample values for six illustrative parameters (βh,βm,ρ,α,andθh) chosen via Latin Hypercube Sampling.
(a)
(b)
(c)
(d)
(e)
(f)


## 4. Discussion

A mathematical model for the transmission dynamics of lymphatic filariasis with treatment for those displaying elephantiasis symptoms is presented as a system of differential equations. The reproduction number and the equilibria are computed and analysed. The model is shown to exhibit backward bifurcation, where a stable disease-free equilibrium coexists with a stable endemic equilibrium when the reproduction number is less than unity. The presence of backward bifurcation makes the classical requirement of having the reproduction number less than unity necessary but insufficient to control the epidemic. Analysis of the reproduction number suggests that treatment of elephantiasis cases has some impact on reducing the spread of lymphatic filariasis infections. However, sensitivity analysis (Monte Carlo simulations) tends to give a better picture about the relationship between the reproduction number and the treatment factor, as they show that the treatment factor reduces the reproduction number but not to levels necessary for disease elimination. This result suggests that effective lymphatic filariasis control requires strategies beyond elephantiasis treatment only. The model presented is not exhaustive; it can be extended to incorporate chemoprophylaxis for the exposed. Given that lymphatic filariasis and malaria are both mosquito-borne infections, a coinfection model of the two would be another worth extension.



---

*Source: 930130-2012-10-17.xml*