A mathematical model of the dynamics of lymphatic filariasis in Caraga Region, the Philippines

Despite being one of the first countries to implement mass drug administration (MDA) for elimination of lymphatic filariasis (LF) in 2001 after a pilot study in 2000, the Philippines is yet to eliminate the disease as a public health problem with 6 out of the 46 endemic provinces still implementing MDA for LF as of 2018. In this work, we propose a mathematical model of the transmission dynamics of LF in the Philippines and a control strategy for its elimination using MDA. Sensitivity analysis using the Latin hypercube sampling and partial rank correlation coefficient methods suggests that the infected human population is most sensitive to the treatment parameters. Using the available LF data in Caraga Region from the Philippine Department of Health, we estimate the treatment rates r1 and r2 using the least-squares parameter estimation technique. Parameter bootstrapping showed small variability in the parameter estimates. Finally, we apply optimal control theory with the objective of minimizing the infected human population and the corresponding implementation cost of MDA, using the treatment coverage γ as the control parameter. Simulation results highlight the importance of maintaining a high MDA coverage per year to effectively minimize the infected population by the year 2030.

Vector control measures such as spraying and use of bednets are generally not a supported strategy for most LF endemic areas in the country due to their underwhelming impact in reducing host-vector contacts. This inefficiency can be attributed to the exophilic, exophagic and day-biting behaviours of Aedes poicilius mosquitoes, which are the principal vectors of LF in the Philippines [32]. These mosquito vectors are known to breed on the leaf axils of abaca and banana plants; hence, any form of vector control that aims to destroy their breeding sites may result to the loss of livelihood of farmers and field workers in these areas [33]. However, in areas wherein Anopheles mosquitoes are the main vectors of LF, vector control measures can enhance the impact of MDA in interrupting the transmission [7].
As of 2018, LF threatens around 893 million people in 49 countries worldwide and a global baseline estimate of 36 million individuals are suffering from chronic disease manifestations [8]. This is a remarkable progress compared to the 1.3 billion people at risk and 120 million people infected with 40 million disfigured and incapacitated in 2000 [34]. However, the magnitude of these numbers reflect how far we are still from our goal of eliminating LF as a global public health problem. Because of this, the WHO has reset its target for elimination of LF to the year 2030 from its initial target of 2020. According to the new NTD roadmap [35], at least 58 out of 72 (81%) endemic countries need to be validated for the elimination of LF as a public health problem by 2030. Unfortunately, the Philippines is one of the countries that are yet to eliminate LF as a public health problem with approximately 5.3 million Filipinos in 6 endemic provinces still requiring MDA as of 2018, which is the DOH's national elimination target date for LF [1,36].
Three general simulation models of LF transmission and control (LYMFASIM [37], EPIFIL [38,39] and TRANSFIL [40]) are currently being used to support policy making and designing elimination strategies for LF. LYMFASIM is an individual-based modelling framework wherein different models for the transmission and control of LF can be composed by choosing different assumptions and by varying parameter values [37]. EPIFIL is a coupled partial differential equation and ordinary differential equation model describing patterns of LF infection and associated diseases through the changes in the human and mosquito populations over age and time [38]. TRANSFIL is described to be the individual-based stochastic equivalent of EPIFIL [40]. All three models account for the human, mosquito and parasite dynamics in the transmission and control of LF. Some studies on these simulation models are discussed in [41][42][43][44][45][46].
Several mathematical models of LF have also looked into the dynamics of the disease by considering the interaction of the human and mosquito populations. In 2009, Supriatna and Soewono investigated the long-term effects of targeted medical treatment in the disease dynamics in Jati Sampurna, West Java, Indonesia using an SI-SI model [47]. Supriatna extended this model by considering an additional human compartment consisting of latent individuals to account for the 'delay' in the infection period [48]. Numerical simulations using an SEI-SI model developed by Bhunu and Mushayabasa suggest that treatment of both exposed and infected humans might lead to a more effective control of the disease compared to the treatment of infected population alone [49]. This model is extended by Bhunu to assess the potential of pre-exposure vaccination and chemoprophylaxis in the control of the disease [50]. In a 2016 study, Iddi et al. assessed the impact of MDA, health education campaigns and the vector control sterile insect technique on the transmission dynamics of LF using a deterministic model [51]. In 2017, Mwamtobe et al. proposed an SEI-SI LF model with quarantine and treatment as control strategies [52]. In these studies, the parameter values are either assumed or obtained from existing malaria models.
Herein, we develop a mathematical model of LF to investigate how mass drug administration impacts the disease dynamics in the Philippines. Compared to the existing models which mostly explored scenarios wherein the treatment is given only to the infectious population, we consider a more realistic representation of MDA wherein the antifilarial drugs are given to all eligible individuals in the population, infected and uninfected alike. Hence, the treatment affects not only the dynamics of the infectious human population but also those who are infected but not yet infectious. Since this study is Philippine-specific, we also use Philippine filariasis data to estimate model parameters to have more meaningful insights. The model considers no intervention on the vector population due to the reasons mentioned previously.
Ultimately, we aim to use optimal control theory (OCT) to determine the best implementation strategy for MDA until 2030 that will lead to a reduction of the infected population with a minimized MDA implementation cost. There are various applications for OCT in many different fields. For instance, in modelling cancer dynamics, OCT can be used to determine the optimal treatment regimen which will minimize tumour density and drug side-effects over a defined period of time [53]. In epidemic modelling, OCT can be used to determine optimal vaccination schedule which will lead to the elimination of the epidemic in the population [54]. In modelling disease transmission and control, OCT royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 can provide insights on the optimal intervention strategy with the least implementation cost [55]. Other studies with applications of optimal control theory can be found in [56,57] (tuberculosis), [58] (dengue) and [59] (HIV/AIDS). An in-depth discussion on optimal control theory and a number of insightful applications are provided by Lenhart et al. [60].
To the authors' knowledge, this study provides the first mathematical model of LF in the Philippines. Although the DOH begun MDA in 2001, we still see Mf prevalence rates of greater than 1% in some areas in the Philippines. This motivates us to study the transmission dynamics of LF and suggest ways to accelerate the elimination of the disease especially in the remaining endemic areas in the country. This study aims to assist the DOH and similar programmes in the region in designing more effective and cost-efficient implementation approaches for MDA fit for each endemic area, to achieve LF elimination in the whole country in the near future.
The rest of the paper is organized as follows. In the Methods section, the proposed model for LF transmission is described. Information on the epidemiological data and the parameter values and initial conditions used are also presented in this section. In the Results and discussion section, stability analysis of the steady state solutions of the model is presented. Results of the sensitivity analysis using the Latin hypercube sampling and partial rank correlation coefficient method are also discussed. The obtained results provide information on the critical model parameters with respect to the infected population, which guide the parameter estimation using the available filariasis data from the Philippine DOH. In this section, numerical simulations on the application of optimal control theory using the forward-backward sweep method are also presented. We investigate the efforts in minimizing the number of infected individuals and the corresponding implementation cost of MDA. We summarize the results of our work and recommend possible extensions of this research in the last section.

Mathematical model of lymphatic filariasis
Filarial parasites that cause LF need two host species to complete their five-stage life cycle: a definitive host (humans or animals) wherein the development from third-stage larva (L3) to adult worm and the reproduction of microfilariae occurs, and an intermediate host (mosquito) wherein the development from microfilaria to L3 occurs. The mosquito also acts as a vector of the parasite that physically carries and transmits infective larvae from one human to another. Without one host or the other, disease transmission will not be sustained in the population. Mosquitoes are infectious to humans if they harbour third-stage larval parasites and humans are infectious to mosquitoes if they harbour microfilariae [61].
The LF transmission dynamics is summarized as follows. When an infected mosquito bites an uninfected human, the infective L3 are introduced onto the skin. The infection in humans begins when these parasites enter the human body through the mosquito bitewound and migrate to the lymphatics where they develop to maturity within 6-12 months [62]. Adult worms reproduce and fecund female worms release thousands of microfilariae. These microfilariae travel through the lymphatic channels into the blood stream where they are taken up by a mosquito through a blood meal. Within 10-12 days, the ingested microfilariae move from the mosquito's gut to its thoracic cavity where they mature to infective L3. At this point, the L3 larvae migrate to the mosquito's proboscis, and the transmission cycle continues when the infected mosquito bites an uninfected human [61].
We propose a deterministic model of LF transmission involving the interaction of the human and mosquito populations. First, we assume that in an LF endemic area, the human population can be categorized into three epidemiological classes based on each individual's infection level: uninfected U h (t), latent L h (t), and infectious I h (t). The latent stage accounts for the development of infective L3 larvae into fecund adult worms. Hence, individuals in the latent stage are considered infected but not yet infectious. Meanwhile, all infected individuals who are able to transmit the infection are in the infectious class. Thus, the total human population is represented by N h (t) = U h (t) + L h (t) + I h (t). On the other hand, mosquitoes can only be either uninfected U v (t), or infected I v (t). Here, the latent stage is considered negligible; thus, infected vectors are assumed to be infectious. We note that the model considers only the female mosquito population since only adult female mosquitoes contribute to the transmission. Hence, the total mosquito population at time t is defined as N v (t) = U v (t) + I v (t). Moreover, the model considers one species of worm and one species of mosquito.
We assume that recruitment into both human and mosquito populations is only by birth. Moreover, since there is no vertical transmission of the infection (i.e. infected pregnant mothers cannot pass the infection to royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 their offspring), all new members of both populations enter their respective uninfected classes at per capita rates b h and b v . Since there is no disease-induced death, both humans and mosquitoes leave their respective population through natural death at respective per capita rates δ h and δ v .
The model also assumes homogeneous mixing of the human and mosquito populations. That is, each uninfected individual in the population has an equal probability of being bitten by an infected mosquito and each uninfected mosquito has the same probability of biting an infected human. Now, if β is the average number of bites per mosquito per unit time, then there are β(N v (t)/N h (t)) mosquito bites per human per unit time. However, only a fraction of these bites, which come from infected mosquitoes, are potentially infective to humans. Further, only a proportion of the potentially infective bites actually result to a successful transmission. Thus, we define the force of infection from mosquitoes to humans, λ vh (t), as the product of the average mosquito bites per human per time β(N v (t)/N h (t)), the probability that the mosquito is infected I v (t)/N v (t), and the probability of successful transmission of infection θ vh , i.e.
Without treatment, the force of infection from humans to mosquitoes, defined as is just the product of the average number of bites per mosquito per unit time β, the probability that the human is infectious I h (t)/N h (t), and the probability of successful transmission of infection θ hv . It is known that the antifilarial drugs given during MDA can instantaneously kill the microfilariae in humans and can halt the reproduction of adult worms, thus reducing the probability of transmission for a significant amount of time [63]. So, if we define p as the proportion of reduction in transmission due to treatment, then the transmission probability becomes θ hv (1 − p). Hence, with treatment, the force of infection is given by Now, suppose only a proportion γ of the human population is given treatment. Then, the total force of infection from humans to mosquitoes is computed using equations (2.1) and (2.2) as follows: Thus, upon sufficient bites to infectious humans, uninfected mosquitoes move to the infected class at the rate λ hv (t). Similarly, after adequate infective mosquito bites, uninfected humans move to the latent class at the rate λ vh (t). Note that the force of infection from mosquitoes to humans remains the same with or without treatment since the antifilarial drugs ingested by humans do not directly affect the mosquitoes. Owing to the growth and development of parasites in the human body, individuals in the latent class will eventually become infectious at the rate α. Since there is no permanent nor temporary immunity, the latent and infectious humans given treatment move back to the uninfected class at the rates r 1 and r 2 , respectively. We note that r 1 ≥ r 2 , since the level of infection (i.e. density and developmental stage of parasites in the body) of the latent individuals is assumed to be lower than those who are infectious. We assume that the parameters b h , b v , δ h , δ v , θ vh , θ hv , β, α are strictly positive real numbers. We also assume that the antifilarial drugs are effective; thus, r 1 , r 2 are strictly positive, and p is within (0, 1]. Since the treatment coverage γ is a controllable parameter, we assume that it can take any value within [0, 1], where γ = 0 implies that no one in the population is given treatment while γ = 1 implies that the whole population is given treatment. We also assume that all individuals given treatment actually ingest the antifilarial drugs. Figure 1 provides a graphical interpretation of our model. Based on our model descriptions and assumptions, the proposed model is governed by the following system of differential equations: (2:3) royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 where

Epidemiological data
The available filariasis data from the Philippine Department of Health's Field Health Service Information System (FHSIS) Annual Reports [64] were compiled for each Philippine province and region covering the years 2009-2018. The 2009-2018 data report annual prevalence rates, computed as the percentage of microfilariae-positive individuals over the total cases examined in the area. According to the DOH FHSIS reports, 5 (out of 5) provinces and 3 (out of 6) cities in the administrative region of Caraga have been endemic for LF [64]. As of 2018, one province (Surigao del Norte) and one city (Surigao City) in Caraga Region are still implementing MDA for LF while the other four provinces have already been declared LF-free (Agusan del Sur and Dinagat Islands in 2010, Surigao del Sur in 2013, and Agusan del Norte in 2015). The 2018 DOH FHSIS report also suggests that there are still confirmed cases of LF in the three provinces and 1 city previously declared LF-free [64]. Hence, the Caraga filariasis data are used in our simulations. From the recorded total human population and annual prevalence rates, the infectious human population per year was estimated as the proportion of the total population that are microfilariae positive (i.e. I h = total population × prevalence rate). Table 1 gives a summary of our dataset.

Parameter values and initial conditions
The human natural death rate δ h , measured in weeks, is computed as the inverse of life expectancy, which is calculated by subtracting the median age of the Caraga population in 2010 [65] from the recorded life expectancy at birth in Caraga in 2010 [66]. An approximate for the obtained value is δ h = 0.00042 ( per week). Using the obtained value and the Caraga population data for the years 2009-2018 in table 1, we estimate the human birth rate b h from the exact solution of Solving this numerically using the built-in Matlab function lsqcurvefit, we find that the birth rate is approximately b h = 0.0006 ( per week). Meanwhile, since LF and dengue are transmitted by mosquitoes of the same genus (Aedes), the parameter values for the mosquito death rate δ v and mosquito biting rate β are obtained from a dengue model presented by de los Reyes V & Escaner IV [67]. Here, we assume that The progression from L h to I h (1/α) is usually between 6 and 12 months [62], but we fix α = 0.0288 which is about eight months, following a study by Jambulingam et al. [44]. In [39], Norman et al. set a value of 1.13 × 10 −4 for the proportion of L3 filarial parasites entering a host which develop into adult worms. In our model, this value is assigned to the transmission probability from vector to human, θ vh . In the same paper, the proportion of mosquitoes which pick up infection when biting an infected host was assigned a value of 0.37, which we set to be the transmission probability from human to vector, θ hv .
The treatment coverage γ = 0.619 is computed as the average of the recorded MDA programme coverages in Caraga from 2010 to 2018 obtained from the DOH FHSIS data [64]. The treatment rates r 1 and r 2 in table 1 are defined as the inverse of the average duration of the treatment of humans in compartments L h and I h , royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 respectively, i.e. the time it takes for an infected person to be treated and to move back to the uninfected class. These parameters are estimated using the available filariasis data in Caraga.
The parameter p represents the proportion of reduction in transmission due to treatment. From what has been inferred from the literature, this parameter depends on many factors such as the percentage of microfilariae and adult worms killed due to treatment, the reduction in Mf production, and the duration of the sustained reduction in the density of these parasites in the human body following treatment. Here, a value of 0.6 is assigned to the parameter p. We show in the next section that the infected human population L h + I h is not sensitive to the parameter p (figure 2).
In our numerical simulations, we use the 2009 population data from Caraga as initial conditions. Since only the total human population N h (0) and the infectious population I h (0) are available, it is assumed that the remaining human population, N h (0) − I h (0), can be apportioned between the uninfected U h (0) and latent L h (0) populations using a partitioning parameter m ∈ [0, 1] in the relation  where the obtained values are rounded off to the nearest whole number to be biologically consistent with the human population. This parameter m is also estimated along with r 1 and r 2 using the available filariasis data from Caraga. Meanwhile, the initial values for the mosquito population are assumed to be U v (0) = 1 000 000 and I v (0) = 100 000. Table 2 lists the parameter values of the LF transmission model.

Mathematical analysis of the model
Since all model parameters and state variables are assumed to be non-negative, the total populations are also non-negative, i.e.
The model system in equation (2.3) is mathematically and epidemiologically well posed in C.
3) has two steady-state solutions: , and : Since E 0 is a steady-state solution wherein the infected populations are zero, E 0 is referred to as the disease-free equilibrium (DFE) solution. The steady-state solution E Ã is known as the endemic equilibrium solution since the infection is constantly maintained in the population. The derivations of E 0 and E Ã are discussed in detail in appendix A.
The basic reproduction number, R 0 , is a threshold parameter that is used to assess whether or not a disease will invade a population. In theory, if R 0 , 1, then each infected person produces less than one new infected individual in their entire period of infectiousness, which implies that the infection will not be sustained in the population. On the contrary, if R 0 . 1, then each infected individual infects more than one person implying that the disease will eventually invade the population. Here, the basic reproduction number is computed using the method of next-generation matrix, a method introduced by Diekmann et al. [68]. The value of R 0 for the model system (2.3) is presented in the next theorem. For the proof, we refer readers to appendices B and C.
Theorem 3.2. The basic reproduction number for the model system (2.3) is Moreover, the DFE E 0 is locally asymptotically stable when R 0 , 1 and unstable when R 0 . 1. The endemic equilibrium E Ã is locally asymptotically stable when R 0 . 1.

Sensitivity analysis and parameter estimation
Sensitivity analysis is a tool used to identify and rank critical input parameters in a model with respect to their impact on the reference model output. An input parameter is said to be influential if small variations of its value result to significant changes in the model output. Based on the result of sensitivity analysis, one can have an idea which input parameters need to be assigned accurate values (i.e. the most influential parameters) and which ones can be roughly estimated (i.e. the less influential parameters). In this study, we used one of the most efficient and reliable global sensitivity analysis techniques-the partial rank correlation coefficient (PRCC) method combined with the Latin hypercube sampling (LHS) technique for the sampling of parameters [69,70]. Numerical simulations for sensitivity analysis using the LHS/PRCC method were carried out using modified versions of the PRCC Matlab codes presented by Marino et al. [69]. The sign and magnitude of the PRCC values, which lie between −1 and +1, characterize the qualitative relationship between the model input and the model output. A positive PRCC implies a positive correlation between the input and the output; that is, an increase in the model input will result in an increase in the model output. On the other hand, a negative PRCC implies a negative correlation wherein an increase in the model input causes a decrease in the model output, and vice versa. The absolute value of a PRCC measures the importance of the model input to the relative model output. The greater the magnitude of a PRCC value, the greater the impact of the input to the output. The p-value of the PRCC indicates the statistical significance of the value. In most cases, a parameter is said to be sensitive or influential to the model output if the magnitude of the PRCC is greater than or equal to 0.5 and the corresponding p-value is less than 0.05 [69]. We used these criteria to determine the influential parameters in our model.
We tested the sensitivity of the model parameters β, α, γ, p, r 1 , r 2 , m and a dummy parameter, with respect to the infected human population L h + I h . A dummy parameter is included to test the robustness of the model as in [69]. That is, since it does not appear in the model equations and does not affect the model in any other way, the dummy parameter should be assigned a sensitivity index of zero in the simulations. The sample size is set to N = 10 000 and the time points of interest are every year for a duration of 10 years or every after 52 weeks for T = 520 weeks.
We fix the values of b h and δ h based on the calculations using the available data. We also fix the values of b v , δ v , and the transmission probabilities θ vh and θ hv , following the results in [39,67].
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 Since there is limited information on the values for α, γ and p, a uniform distribution is assigned to their corresponding parameter ranges which we set to +50% of their respective baseline values given in table 2. On the other hand, since β has already been estimated in [67], we used a normal distribution with the nominal value in table 2 as the mean and a standard deviation of 0.1. The range of values for the treatment rates r 1 and r 2 is from 0.0019 to 0.5, or from two weeks to 10 years, while the range of values for the partitioning parameter m is from 0 to 1. Both ranges are assigned a uniform distribution.
The result of our simulations is shown in figure 2, where each coloured bar graph corresponds to a particular time point. Based on our criteria for determining critical parameters, the results suggest that the parameters γ, r 1 and r 2 are the most influential parameters throughout the running time of 10 years or 520 weeks. It can be observed that γ, r 1 and r 2 have a negative correlation to the model output L h + I h . This implies that an increase in the treatment coverage and faster treatment rates will most likely result in a decrease in the total infected human population. We also observed that all three parameters have high magnitude of PRCC values for the first five years, after which the magnitude of the PRCC values of γ and r 1 decreased below 0.5. The PRCC values of r 2 decreased below 0.5 in the last two years.
In our parameter estimation, we opted not to include the treatment coverage γ since there are already recorded values for this in the DOH FHSIS filariasis data [64]. This left us with the parameters r 1 and r 2 , which we estimate together with the partitioning parameter m, using the available population data from Caraga Region [64] covering the years 2009-2018.
From the DOH FHSIS filariasis data in Caraga Region in table 1, we have 10 data points for I h (2009-2018) with the 2009 data as the initial condition. The parameters r 1 , r 2 and m are estimated by minimizing the mean of the squared difference between the available data and the model output at the corresponding time point using the Matlab function lsqcurvefit.
To establish unbiased choice of initial conditions for the parameter estimation, the LHS method is used to generate 1000 sets of initial guesses for the estimates. For each set of initial guesses, parameters are estimated using lsqcurvefit. The 'best' estimate is then computed as the average of the 1000 estimates from 1000 distinct initial guesses. The simulation results, along with the lower bound and upper bound set for each of the three parameters for the Matlab simulations are given in table 3.
Using the estimated value for the partitioning parameter m = 0.853284 and the 2009 filariasis data from Caraga in table 1, we can now solve for the initial populations for U h and L h in equation (2.4). Meanwhile, the estimated value for the treatment rate r 1 is 0.430848 which means that the treatment period 1/r 1 is about 16 days, and the estimated value for the treatment rate r 2 is 0.010038 which implies that 1/r 2 is about 2 years. Recall that the progression from the latent stage to the infectious stage lasts for about 6-12 months [62]. Hence, the obtained value for r 1 highlights the effectiveness of MDA in preventing new infections. The obtained value for r 2 explains that although transmission of L3 larvae from mosquitoes to humans is inefficient, the treatment of infectious individuals takes a long time because of the 5-year reproductive life span of adult filarial worms [61]. This highlights the importance of maintaining a high treatment coverage since majority of the infected population are asymptomatic. Figure 3 depicts the available filariasis data plotted with the corresponding model fit from parameter estimation.
Parameter bootstrapping is also used to investigate the reliability of our estimated values. Bootstrapping is a statistical technique used to quantify uncertainty and construct confidence intervals in parameter estimates [71]. Following the algorithm presented by Chowell [71], 1000 samples of synthetic datasets are generated from the best-fit model by assuming a Poisson error structure. Parameters are then re-estimated from each synthetic dataset using lsqcurvefit to obtain a new set Table 3. Results of parameter estimation using lsqcurvefit with the bounds used. The estimate is computed as the average of the estimates from 1000 distinct initial guesses generated by the LHS method. of parameter estimates. Numerical simulations for parameter bootstrapping were carried out using modified versions of the sample Matlab codes presented by Chowell [71]. Figure 4 shows the parameter bootstrapping results with the means, standard deviations and 95% confidence intervals for each parameter. We observe low standard deviations and narrow ranges of the 95% confidence intervals for all the parameters. Moreover, the estimates obtained for m, r 1 and r 2 (table 3) all lie within their corresponding confidence intervals.

Optimal control strategy
In implementing MDA programmes, the treatment coverage can vary through time depending on the efforts of the government. In the LF transmission model, this translates to the constant parameter γ becoming a function of time u(t). Since the main purpose of MDA is to limit and reduce the infected humans by interrupting the transmission of infection, the proposed optimal control problem is to minimize the number of infected (both latent and infectious) hosts as well as royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 the corresponding cost of implementing MDA. Here, minimizing the infected mosquito population is not included in the objective since there is no control directly affecting the vector population. From our simulations, we aim to gain insights about the optimal MDA coverage over a duration of T years, T = (t f − t 0 ) > 0, such that the objective is satisfied. Hence, our objective functional to be minimized is where t 0 and t f are taken as 2018 and 2030, respectively. Here, c is a weighting parameter associated with the implementation cost of MDA. We note that the parameter c does not represent per se the actual monetary cost of implementing MDA. Instead, it is a constant parameter that balances the size and importance of each term in the integrand. That is, if c is too high, more importance will be given in reducing the cost of implementation of MDA compared to the infected population. On the other hand, if c is low, the minimization problem will put equal importance in minimizing the infectious population and the cost of implementation. It is assumed that the control is a quadratic function to represent nonlinear implementation costs. Our goal is to find optimal u Ã satisfying with the initial conditions U h (0) = U h,0 , L h (0) = L h,0 , I h (0) = I h,0 , U v (0) = U v,0 , I v (0) = I v,0 , and such that u min ≤ u(t) ≤ u max . Pontryagin's maximum principle [72] transforms the optimal control problem into a problem that minimizes a Hamiltonian H pointwise with respect to the control u(t). Using and λ(t) = [λ 1 , …, λ 5 ](t), the Hamiltonian H is formed as Applying Pontryagin's maximum principle, we obtain the following result. The proof of this theorem can be found in appendix D.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 Theorem 3.3. There exists optimal control u Ã (t) minimizing the objective functional J(u(t)) over V ¼ {u min u u max ; u [ L 2 ð2018; 2030Þ}. Given this optimal solution, there exist adjoint variables λ 1 (t), …, λ 5 (t), which satisfy with transversality conditions λ i (t f ) = 0, for i = 1, …, 5. Furthermore, the optimality equation is as follows: Since our compiled population data in table 1 suggest that there are still infectious individuals in Caraga in 2018, we look at the disease dynamics after 2018 to obtain insights on the optimal strategy for implementing MDA to eliminate LF in the population by 2030. Hence, the choice of t 0 = 2018 and t f = 2030 in the objective functional.
Numerical simulations were carried out in Matlab using the forward-backward sweep method [60]. We set u min = 0 to account for the possibility that MDA is stopped after 2018, and u max = 0.95 as we recognize that a 100% treatment coverage is quite difficult to achieve in reality. We also considered different values for the weighting parameter c to investigate how the variations in the MDA implementation cost affect the optimal control solution. The obtained optimal control solutions for c = 10 i , i = 0, 2, 4, 6, were compared with the no control solution (u(t) ¼ 0, 8t) and the constant control solution (u(t) ¼ 0:619, 8t). Figures 5 and 6 show the obtained control profile and the corresponding effect in the dynamics of the infectious population I h (t).
Observe from figure 6 that if MDA is stopped after 2018, the decrease in the infectious population over a span of 12 years is only about 15%. On the other hand, if the current control is maintained for another 12 years the infectious population can decrease further by 98%.
One apparent conclusion from the optimal control solutions in figure 5 is that the MDA coverage is inversely proportional to the implementation cost. That is, the lower the cost, the higher the MDA coverage. Consequently, a higher MDA coverage leads to a lower number of infected individuals at the end time as illustrated in figure 6. Here, we compare the optimal control results for different values of c to the no MDA strategy. As recorded in table 4, with c values equal to 10 0 , 10 2 or 10 4 , the I h population can be reduced by at least 91% at the end of the control period. When c = 10 6 , we observe that the optimal control problem gives more emphasis on minimizing the MDA implementation cost compared

Conclusion
LF remains a public health problem in the Philippines with millions of Filipinos still at risk of being infected as of 2018. In this work, we developed a mathematical model of LF transmission involving the interaction of human and mosquito populations. Using this model, we investigated how the implementation of the annual MDA for LF affects the disease dynamics in the LF-endemic region of Caraga. Sensitivity analysis using the LHS/PRCC method showed that the infected human population is most sensitive to the treatment coverage (i.e. how much of the population receives treatment) and the treatment rates (i.e. how effective the antifilarial drugs are in reducing the parasite density in infected humans). This highlights the importance of strategic MDA implementation, and efficient data gathering and data reporting to obtain more accurate simulation results and produce more realistic and relevant insights. The treatment parameters were estimated using the available filariasis data in Caraga obtained from the Philippine DOH's FHSIS. The estimated values emphasize the importance of MDA in preventing new infections.
Parameter bootstrapping showed small variability in all the parameter estimates, indicated by the low standard deviations and narrow confidence intervals. The application of optimal control theory highlighted the importance of maintaining a high treatment coverage per year to effectively minimize the infected population by the year 2030. To achieve this, it is  Figure 6. The corresponding effect of the controls in the dynamics of the infectious population I h . The y-axis is in log scale to better illustrate the difference in the dynamics towards the end time for different values of c. Table 4. Per cent reduction in the infectious population I h when the optimal control using implementation cost c is applied. Per cent reduction is the relative difference in percentage between the number of infectious population in 2030 when optimal control is applied and without MDA. important that there is an effective and systematic implementation of MDA from the national level down to the implementation units. Activities such as the distribution of antifilarial drugs, collection of data, and diagnostic examinations must be enhanced. To address the problems with MDA acceptance [73], educational activities must be conducted before the distribution of antifilarial drugs to inform the people and help them understand how participating in MDA will benefit them and their communities. Field workers and volunteers must also be well informed about the MDA programme and the disease. Most importantly, since LF is closely related to poverty, eliminating the disease requires the active commitment of the government to improve the lives of its poorest and most vulnerable citizens. This does not only include the distribution of treatment against LF and other NTDs but also the provision of better living conditions by giving people better opportunities for employment, proper education, clean water and other basic human needs. One limitation of our study is the uncertainty in the values of our model inputs, both initial conditions and model parameters. Since the model is applied specifically to Caraga, the parameter values must be accurate to well represent the transmission dynamics of LF in the region. However, majority of the existing studies on LF in the Philippines are outdated and may no longer give accurate information on the disease dynamics at present. We urge the Philippine DOH to have a systematic, correct, and timely data gathering that may be relevant to the study of not only LF but also other diseases. We also encourage local researchers to look into LF in the country. One possible topic to explore is the dynamics of Aedes poicilius mosquitoes, or other known mosquito species that transmit LF in the Philippines, and review values of relevant parameters. Both clinical and field studies on the efficacy of the antifilarial drugs, or a comprehensive study on the general effects of MDA in reducing transmission in a population under varying circumstances, could also be helpful for this study. A possible extension of our model is the inclusion of the parasite dynamics, such as parasite aggregation, to better capture the transmission processes of LF. Future works could also focus on calibrating the model to incorporate heterogeneity in susceptibility and infectivity of individuals in the population. We also recommend cost analysis of MDA in the remaining endemic areas, such as the work of Amarillo et al. in Sorsogon, Philippines in 2009 [74].
Data accessibility. The filariasis data of the total population, total cases examined, cases found positive, and prevalence rate for Caraga Region, the Philippines, presented in table 1 are obtained from the Philippine Department of Health's Field Health Service Information System Annual Reports [64].
Authors' contributions. P.K.N.S. participated in the conceptualization, methodology, and formal analysis of the study, curated the data, carried out the simulations and computations, drafted the manuscript, acquired funding, and participated in the critical review of the manuscript; V.M.P.M. and R.G.M. participated in the conceptualization, methodology, and formal analysis of the study, drafted the manuscript, acquired funding, supervised the study, and participated in the critical review of the manuscript; V.Y.B.J. participated in the conceptualization and formal analysis of the study, supervised the study, and participated in the critical review of the manuscript. All authors participated in the final review and editing of the manuscript, and gave final approval for publication and agree to be held accountable for the work performed therein.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 where we look for the corresponding solution to We have the following: We first solve the state variables in terms of λ vh , then express λ vh in terms of the model parameters.

Using equation (A 4) and equation (A 5), we get
: Moreover, if we plug-in equation (A 2) to equation (A 3), we get : Plugging-in equation (A 7) and equation (A 2) to equation (A 1) and doing a few algebraic manipulations, we have : Using equation (A 8) in equation (A 2), we get : Similarly, using equations (A 9) and (A 2), we obtain : In the same manner, we find an expression for λ hv using equation (A 10): Then we have and : royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 Using equation (A 12) and equation (A 4), we get : Furthermore, plugging-in equation (A 13) to equation (A 6), we get : Plugging-in equation (A 12) to equation (A 11) and after doing a few algebraic computations, we obtain and Here, we have two possible cases: either λ vh = 0 or Bλ vh − A = 0.
Case 2. Suppose Bλ vh − A = 0 and λ vh ≠ 0. Then : We denote the above as l Ã vh . We have the solution , and : Note that E Ã is biologically feasible if l Ã vh . 0. Our assumptions on the positivity of parameter values imply that B is always positive. It can also be shown that A > 0, i.e.

Appendix B. Computation of the basic reproduction number
Following van Driessche & Watmough [75], we first rearrange our system such that the first n entries contain the infected populations. We have _ royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 Since there are three compartments for the infected population in our model, we have n = 3. Now, _ X can be expressed as where F i (X) is the rate of appearance of new infections in compartment i, and V i (X) ¼ V À i (X) À V þ i (X) is the rate of transfer of individuals into compartment i by all other means (V À i (X)) and out of compartment i (V þ i (X)), for i = 1, …, 5. In our model, we have : From here, we obtain the matrices where X 0 is a disease-free equilibrium and i, j ∈ {1, 2, 3}. Diekmann et al. [68] defined R 0 as the spectral radius, or the dominant eigenvalue, of the matrix FV −1 . Following this, we have and evaluated at the DFE, We also obtain Now, solving for the product of F(E 0 ) and V −1 (E 0 ), we have Taking the largest eigenvalue of the above, we have s :

Appendix C. Stability analysis of the steady-state solutions
To study the stability of our equilibrium solutions, we first linearize the system (2.3) near the steady states and calculate the corresponding Jacobian.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 Àd v   2   6  6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  7  5 : The corresponding characteristic equation is To prove that E 0 is locally asymptomatically stable, we need to show that the eigenvalues of (C1), i.e. the roots of the characteristic equation above, are all negative. By assumption, δ h , δ v are strictly positive, thus we only need to show that the roots of the cubic polynomial By the Routh-Hurwitz criterion for stability [76], the roots of equation (C 2) are negative if the coefficients of the cubic polynomial above are all positive and a 1 a 2 > a 0 a 3 . It follows from our assumption on parameter values that Note that a 3 > 0 if and only if R 0 < 1. Indeed, . 0 1 À R 2 0 . 0 R 0 , 1: Now, to show that a 1 a 2 > a 0 a 3 , let Then, we can write royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 Thus, a 1 a 2 À a 0 a 3 ¼ ( since A 1 , B 1 , C 1 > 0 and D 1 ≥ 0. Thus, the DFE is locally asymptotically stable when R 0 , 1. On the other hand, when R 0 . 1, the coefficient a 3 < 0. By Descartes' rule of signs [76], this implies that there will be one positive eigenvalue which proves the instability of E 0 . Therefore, the disease-free steady state E 0 is locally asymptotically stable when R 0 , 1 and unstable when R 0 . 1.
To study the stability of E Ã , we look at the eigenvalues of the Jacobian matrix evaluated at E Ã , that is, Àd v : The corresponding characteristic equation is : Since δ h , δ v are strictly positive for all time t, t ≥ 0, we only need to show that the roots of the cubic polynomial are negative. Substituting the values for U Ã h and U Ã v , we obtain royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201965 By the Routh-Hurwitz criterion for stability, the roots of (C 3) are negative if the coefficients of the cubic polynomial above are positive and b 1 b 2 > b 0 b 3 . It follows from our assumption on the parameter values that . 0: Note that b 3 > 0 if and only if . 0: Following a rigorous computation, the above inequality is equivalent to R 2 0 . R 0 which means R 0 . 1 since R 0 . 0. Now, to show that b 1 b 2 > b 0 b 3 , let and : Then, we can write We have . 0 since l Ã vh . 0, A 2 , B 2 , C 2 , D 2 > 0 and E 2 ≥ 0. Hence, the endemic equilibrium E Ã is locally asymptotically stable when R 0 . 1.