Transmission Dynamics of Lymphatic Filariasis: A Mathematical Approach

An epidemiological model for the spread of lymphatic ﬁlariasis, 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 ﬁlariasis 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 ﬁlariasis cases displaying elephantiasis symptoms are put on treatment it will not be able to eradicate the disease. This result suggests that e ﬀ ective control of lymphatic ﬁlariasis may lie in treatment for those displaying symptoms as well as chemoprophylaxis for the exposed.


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 tissuedwelling nematodes of Brugia malayi, Brugia timori, and Wuchereria bancrofti species [1,3] and is estimated to affect about 120 million people worldwide [4][5][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][14][15][16][17][18][19], to mention just a few, but only a few have looked into lymphatic filariasis [20][21][22][23][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.

Model Description
The model subdivides the human population based on lymphatic filariasis status into the following subpopulations: susceptible humans S h (t), infected and not showing signs of elephantiasis (exposed) E h (t), infected and displaying elephantiasis symptoms I h (t). Thus, the total human population is given by (1) The mosquito population is divided into the following subgroups: susceptible mosquitoes S m (t) and infected mosquitoes I m (t), so the total mosquito population is given by 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 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 I m (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 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 E h (t) (infected and not showing signs of filariasis). Individuals in E h (t) progress to the stage of showing filarasis symptoms I h (t) at rate ρ. However, some progress to the I h (t) as a result of reinfection at a rate λ m (t). Individuals in I h (t)-class are treated using Diethycarbamzime (DEC) and albendazole drugs at a rate α to move back into the susceptible class S h (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. Unless stated otherwise, values for the parameters in the simulations are given in Table 1 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.
Due to lack of data to calibrate the model and for parameter estimation, other parameter values are assumed  Humans  Table 1. 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: All parameters and state variables for model system (5) are assumed to be nonnegative to be consistent with human and mosquito populations. All feasible solutions of model system (5) are positive and eventually enter the invariant attracting region It is sufficient to consider solutions in Ω. Existence, uniqueness, and continuation results for system (5) hold in this region, and all solutions starting in Ω remain there for all t ≥ 0. Hence, (5) is mathematically and epidemiologically well posed, and it is sufficient to consider the dynamics of the flow generated by the model system (5) in Ω.  Table 1.

Disease-Free Equilibrium and Stability Analysis.
The disease-free equilibrium of model system (5) is given by Following van den Driessche and Watmough [26], we have 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].

Mosquito infection transmissibility
Human infection transmissibility Figure 5: Partial rank correlation coefficients showing the effect of parameter variations on R E using ranges in Table 1. Parameters with positive PRCCs will increase R E when they are increased, whereas parameters with negative PRCCs will decrease R E when they are increased.
The fact that ∂R E /∂α < 0 suggests that increase in treatment of elephantiasis cases results in a decrease of the lymphatic filariasis infections in the population.

Existence and Stability Analysis of the Endemic Equilibrium. The endemic equilibrium is given by
Substituting (10) into the force of infection λ * m in (4), we obtain Also, substituting (10) into the force of infection λ * h in (3), we obtain 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 Clearly A is positive throughout, C > 0 (C < 0) only when R E < 1(R E > 1) and B can be positive or negative depending The analysis of (14) is Theorem 2.
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 S m = x 1 , I m = x 2 , so that N m = 2 n=1 x n and S h = x 3 , E h = x 4 , I h = x 5 , so that N h = 3 n=1 x 2+n . We now use the vector notation X = (x 1 , x 2 , x 3 , x 4 , x 5 ) T . Then, model system (5) can be written in the form dX/dt

ISRN Biomathematics
The Jacobian matrix of system (15) at E 0 is given by from which it can be shown that the reproduction number is 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 R E = 1 and solve for β h , we obtain Note that the linearized system of the transformed equation (15) with the bifurcation point β * has a simple zero eigenvalue. Hence, the Centre Manifold Theory [27] can be used to analyze the dynamics of (15) near β h = β * . It can be shown that the Jacobian of (15) at β h = β * has a right eigenvector associated with the zero eigenvalue given by u = [u 1 , u 2 , u 3 , u 4 , u 5 ] T , where The left eigenvector of J(E 0 ) associated with the zero eigenvalue at β h = β * is given by
, a > 0, then model system (5) undergoes a backward bifurcation at R E = 1, otherwise a < 0 and a unique endemic equilibrium E * guaranteed by Theorem 4.1 [28] is locally asymptotically stable for R 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 R 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.

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 (5) using model parameters in Table 1. Figure 3 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.
In Figure 4, 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.

Sensitivity Analysis.
In order to investigate the effects of variations in R 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. Figure 5 illustrates the PRCCs using R 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 R 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 6 illustrates the effect that varying six sample parameters will have on R E . If the mosquito death rates are sufficiently high, then the disease prevalence can be reduced (R E → 1), but not sufficiently enough to control it (R E < 1). Similarly, high treatment levels result in R E → 1 but are not sufficient enough to have R E < 1 necessary for disease control. Increases in the average number of mosquito bites (β m , β h ) increase R 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.

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 diseasefree 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.