Delays in lymphatic filariasis elimination programmes due to COVID-19, and possible mitigation strategies

Abstract Background In view of the current global coronavirus disease 2019 pandemic, mass drug administration interventions for neglected tropical diseases, including lymphatic filariasis (LF), have been halted. We used mathematical modelling to estimate the impact of delaying or cancelling treatment rounds and explore possible mitigation strategies. Methods We used three established LF transmission models to simulate infection trends in settings with annual treatment rounds and programme delays in 2020 of 6, 12, 18 or 24 months. We then evaluated the impact of various mitigation strategies upon resuming activities. Results The delay in achieving the elimination goals is on average similar to the number of years the treatment rounds are missed. Enhanced interventions implemented for as little as 1 y can allow catch-up on the progress lost and, if maintained throughout the programme, can lead to acceleration of up to 3 y. Conclusions In general, a short delay in the programme does not cause a major delay in achieving the goals. Impact is strongest in high-endemicity areas. Mitigation strategies such as biannual treatment or increased coverage are key to minimizing the impact of the disruption once the programme resumes and lead to potential acceleration should these enhanced strategies be maintained.

In all three timings above, the round due in 2020 is cancelled and the programmes resume one year later. Following the methodology in the main text, we extracted 10,000 simulations from the TRANSFIL model, uniformly, in the prevalence range between 1% and 40%. To summarize results over the prevalence range, we calculated the moving average with a window size of 2,000 and the volatility as the unweighted standard deviation in the same window size.
The results shown in Figure S1 suggest that the timing of the disruption does not affect the average delay of the programme. Note that for disruptions at the start of the programme, cancelling the first round, essentially translates to a one year delay. The uncertainty is however much wider in that scenario, as treatment rounds reduce the variation in mf prevalence across the simulations. Figure S1. Moving average of the delay in reaching 1% mf prevalence (in years) for a wide range of baseline prevalence at the start of the current MDA programme. The red, orange and lime lines illustrate respectively a disruption mid-programme, at the end of the programme and at the beginning of the programme. Shaded colour areas illustrate the standard deviation (volatility).

Sensitivity to vector control use
We conducted additional simulations to evaluate the impact of vector control in the context of the disruptions explored in the main manuscript. We used the model EPIFIL to simulate populations in a setting with Culex transmission treated with DA and assumed that a bednet and an indoor residual spraying (IRS) coverage of 30% was achieved, and was maintained without disruption (even when MDAs are delayed). The implementation of vector control in the model is detailed below in the EPIFIL model description section. We explored two delays in MDA due to the disruption caused by COVID-19: a 12 month delay and a 24 month delay in the delivery of treatment.
Following the methodology in the main text, we extracted 10,000 simulations from the EPIFIL model, uniformly, in the prevalence range between 1% and 40%. To summarize results over the prevalence range, we calculated the moving average with a window size of 2,000 and the volatility as the unweighted standard deviation in the same window size.
The results shown in Figure S2 -left show the EPIFIL results similar to Figure 3 in the main text. Adding low level vector control can potentially reduce the average delays by around 15% and 32% for the 12 and 24 month delay scenarios respectively, Figure S2 -right. It is important to note the baseline in the counterfactual itself (i.e. no disruption scenario) will already change due to the inclusion of vector control, with a decreased prevalence, thus reducing delays by half a year on average (not shown). Figure S2. Moving average of the delay in reaching 1% mf prevalence (in years) for a wide range of baseline prevalence at the start of the current MDA programme in an a Culex-transmitted setting treating with DA, with vector control (right) and without (right). The red and yellow lines illustrate scenarios with a 12 or 24 months programme delay respectively; shaded colour areas illustrate the standard deviation (volatility).

The mathematical model of LF transmission dynamics
We employed a genus specific mosquito-vectored transmission model of LF to carry out the modelling work in this study [1][2][3][4][5][6][7] . Briefly, the state variables of this hybrid coupled partial differential and differential equation model vary over age (a) and/or time (t), representing changes in the pre-patent worm burden per human host adult worm burden per human host the microfilariae (mf) level in the human host modified to reflect infection detection in a 1 mL blood sample the average number of infective L3 larval stages per mosquito (L), and a measure of immunity developed by human hosts against L3 larvae. The state equations comprising this model are: The above equations involve partial derivatives of four state variables (P -pre-patent worm load; W -adult worm load; M -microfilaria intensity; I -immunity to acquiring new infection due to the pre-existing total worm load where W T = W(a,t) + P(a,t)). Given the faster time scale of infection dynamics in the vector compared to the human host, the infective L3-stage larval density in mosquito population is modelled by an ordinary differential equation essentially reflecting the significantly faster time-scale of the infection dynamics in the vector hosts. This allows us to make the simplifying assumption that the density of infective stage larvae in the vector population reaches a dynamic equilibrium (denoted by L * ) rapidly 1,2,5,8,9 . This basic coupled immigration-death structure of the model as well as its recent extensions has been extensively discussed previously 1-3, 5, 8, 9 . The effects of worm patency are captured by considering that at any time t, human individuals of age less than or equal to the pre-patency period, τ, will have no adult worms or mf, and the rate at which pre-patent worms survive to become adult worms in these individuals at a > τ is given by . The term enables us to account for the different establishment and development rates of the incoming L3-stage larvae as adult worms depending on the genus of mosquito vectors as expressed below: for mosquitoes of Anopheline genus, for mosquitoes of the Culicine genus.
In the above, k [ k 0 +k Lin M (a , t ) ] is the shape parameter of the negative binomial distribution on the mf uptake whereas r and are respectively the rate of initial increase and the maximum level of L3 larvae. See Table S2 for the description of all the model parameters and functions.

Mathematical expressions of the functions Parameters
Probability that an individual is of age a π(a) Human age a in month, A0 and B0 estimated from country demographic data 1,5,9 Larvae establishment rate (modified by acquired immunity)

Ω(a,t)
-proportion of L3 leaving mosquito per bite; -the establishment rate 1 -Adult worm mating probability ϕ(W,k) k -negative binomial aggregation parameter 2,5,22 Immunity to larval establishment g1(I) 1 1+cI c -strength of immunity to larval establishment 1,5 Host immunosuppression g2(WT) SC -slope of immunosuppression 1, 5 1 The proportion of L3-stage larvae infecting human hosts that survive to develop into adult worms 2 . 2 The gradient of mf uptake r is a measure of the initial increase in the infective L3 larvae uptake by vector as M increases from 0 2, 9 . 3 The facilitated establishment rate of adult worms due to parasite-induced immunosuppression in a heavily infected human host 4 The initial rate of increase by which the strength of immunosuppression is achieved as W increases from 0 23 . # Note MBR (monthly biting rate) serves as an input to initialize the model, measured as mosquito bites per person per month, the value of which may be obtained from entomological surveys conducted in study sites. In the absence of the observed MBR value, the model has been adapted to estimate it from the community-level mf prevalence data.

Parameter selection and simulation procedure
We employed the Bayesian Melding (BM) procedure to calibrate and estimate the LF models from field data, as outlined in detail in our previous work 2, [5][6][7] . Typically, we begin the procedure by first using the known or assignment of a uniform range for parameter values to generate distributions of parameter priors. We then randomly sample with replacement from these prior distributions to generate 200,000 parameter vectors, which are run using the annual biting rate (ABR) values, if given, for a site to generate model outputs. The model outputs are then melded with age-stratified mf prevalence data by calculating binomial log-likelihoods for each parameter vector. In the resampling step of the BM method, a Sampling-Importance-Resampling (SIR) algorithm is used to perform 500 draws with replacement from the pool of parameter vectors generated as above, with probabilities proportional to their relative log likelihood values. This step selects the parameter vectors which best describe the given mf ageprevalence data. These resampled parameter vectors are then used to generate distributions of variables of interest from the fitted model (eg. age-prevalence curves, worm breakpoints, and infection trajectories following treatments).
Here, we modified our standard Monte-Carlo BM framework for model discovery 2,[5][6][7]24 to provide simulations for the chosen scenarios for two settings (one relevant for most of Africa, with Anopheles-driven transmission and annual treatment with ivermectin and albendazole (IA); and a second representing India-like populations, with Culex-driven transmission and treatment with diethylcarbamazine and albendazole (DA). The explicit aim was to generate at least N = 100,000 parameter vectors which resulted of overall mf prevalence values 1-60% at baseline in 2018. We first randomly sampled n = 100,000 parameter vectors from the assigned uniform parameter priors. We then simulated the endemic equilibrium given each parameter vector and calculated the predicted mf prevalence at baseline. Those parameter vectors whose outputs produced mf prevalence values between 1-60% were accepted while all others are rejected. This sampling and acceptance/rejection procedure was repeated until the total number of accepted parameter vectors (N) was greater than or equal to 100,000. The N posterior parameter vectors were then used to simulate the impacts of MDA and vector control interventions.

Modeling intervention by mass drug administration
Intervention by mass drug administration was modeled based on the assumptions that antifilarial treatment with a combination drug regimen acts by killing certain fractions of the populations of adult worms and microfilariae instantly after the drug administration 25 . These effects are incorporated into the basic model by calculating the population sizes of worms and microfilariae as follows: where dt is a short time period since the ith MDA was administered. During this short time interval, a given proportion of adult worms and microfilariae are instantly removed. The parameters ω and ε are drug killing efficacy rates for the two life stages of the parasite while the parameter C represents the MDA coverage. Apart from instantaneous killing of microfilariae, the drug continues to kill the newly reproduced mf by any surviving adult worms at a rate δ reduc for a period of time, p. We model this effect as follows: We simulated LF intervention by running the model with fixed values of ω, ε, δ reduc , and p for MDA coverage levels given by the scenarios. The first MDA round was implemented in the model by affecting the population sizes of worms and microfilariae from the baseline estimates, and then the intervention is simulated forward in time for a number of years, with subsequent MDA rounds implemented annually or biannually.

Modeling intervention by Vector Control
In addition to MDA, we also modeled the added effect of long lasting insecticidal nets (LLINs) and indoor residual spray (IRS) as described previously 6 . The impact of LLINs with three main actions against mosquito biting was modelled: 1) deterrence from entering the home (efficacy η 1 ), 2) inhibition of their ability to feed on humans (efficacy η 2 ), and 3) killing them (efficacy η 3 ) 26,27 and the impact of IRS with two main actions against mosquito biting was modelled: 1) inhibition of their ability to feed on humans (efficacy η 4 ), and 2) killing them (efficacy η 5 ) 26,27 . To capture these effects, which decay over time as the larvicide efficacy declines exponentially at rate Λ, we adjust the term V/H to be appropriately modified according to the population coverage for LLINs (C LLIN ) and IRS (C IRS ):

LYMFASIM model description and methods
Description of the mathematical model LYMFASIM 1,2 is a stochastic individual-based model for lymphatic filariasis (LF). It is a specific model variant within WORMSIM, a generalized framework for modelling transmission and control of helminth infections in humans 3,4 . LYMFASIM simulates the life histories of individual people and individual worms in a community, and the effects of interventions (e.g. mass drug administration, integrated vector management, bednet use) on transmission and morbidity, while taking into account of the human demography and the complexities of helminth transmission. The model has been described elsewhere and has been applied to support decision making on control and elimination of lymphatic filariasis in different settings 1,2,5-12 .
Mass drug administration (MDA) is simulated by specifying the year and month in which treatment takes place, the efficacy of the applied treatment regimen, the achieved coverage level, and compliance patterns. LYMFASIM assumes that a fraction of the population never participates in MDA (e.g. systematic refusal, related to chronic illness). In addition, LYMFASIM allows the relative compliance to vary between age and sex groups; this mechanism captures transient contra-indications for MDA (e.g. exclusion of young children and pregnant women) and other age-and sex-related behavioral factors driving participation in MDA. Lastly, each individual has a personal inclination to participate in MDA, which is considered as a lifelong property. A stochastic process eventually defines for each individual whether they are treated in a given round, depending on the calculated probability.

Parameter quantification and simulation methods for this study
We previously derived model quantifications for simulating LF transmission by Culex quinquefasciatus in India 2 and by Anopheles species in Africa 8 . These model quantifications differ in the age-structure of the human population and the density dependence in the L3 yield from a blood meal in mosquitoes. For India, we originally developed several model variants 2 : two model variants assumed that people would acquire some form of acquired immunity against infection during their lifetime, as a consequence of being exposed to various stages of infection; one model variant assumed no immunity. The variants with immunity fitted best to the age-specific data that showed a decline in prevalence in older age group. Of these, the 'anti-L3' variant (with immunity triggered by incoming L3 larvae and reducing the probability for these larvae to survive and develop successfully into adult worms) has been used most in later studies [5][6][7] as it seemed most consistent with external antigenaemia prevalence data. As a result of strong acquired immunity, the range of mf prevalences that can be simulated with this model is relatively narrow.
A review of age-prevalence patterns 13 has casted doubt about the assumed role of acquired immunity. In view of these doubts, acquired immunity was not considered to play a role in the Africa model 8 . We now also revert to the model variant without acquired immunity for India. As shown by Subramanian, exclusion of immunity requires some changes in other parameters, notably the success ratio (i.e. the proportion of incoming L3 larvae that survives and develops successfully into an adult worm) and the assumed variation in exposure between individuals 2 . In the anti-L3 immunity model variant, the success ratio is down regulated as a consequence of acquired immunity. In a model without such immunity, we therefore need to assume a lower value. Also, to simulate the relatively low pre-control mf prevalence levels that are common in India without acquired immunity, we need to assume more heterogeneity in exposure with many individuals having a low exposure and a small group of people being highly exposed to mosquito bites. Subramanian proposed a model in which part of the simulated population does not participate in transmission. We take a slightly different approach in this paper, combining relatively low average monthly biting rates (mbr) with high exposure heterogeneity. The monthly biting rate and parameter describing the extent of exposure heterogeneity (i.e. the shape parameter of the Gamma distribution that describes exposure heterogeneity) are varied according to a predefined parameter space as defined below.
Parameter values are listed in Table S3. Assumptions and parameters related to control strategies and treatment efficacy are listed in Table S4. The monthly biting rate (mbr), exposure heterogeneity parameter (k) and external force-of-infection (foi) were varied according to the density plots in Figure S5, in order to generate simulations across a wide range of mf prevalences at baseline, measured in the total population (all ages).

Version and code availability
For this paper, we used WORMSIM version 2.58Ap27, available from https://gitlab.com/erasmusmc-public-health/wormsim.previous.versions/-/blob/master/ wormsim-2.58Ap27.zip. Figure S5. Density plots illustrating the parameter space areas for the monthly biting rate (mbr), parameter describing exposure heterogeneity in the population (k) and external force of infection (foi), as used to generate simulations with baseline mf prevalence in the range from 0 to 40%, for the Africa and India model variants. Proportion removed when maximum population size is reached 5% 5% N/A

Anopheles model Culex model Source / remarks
External force-of-infection at start of burn-in period 2 2 N/A Duration of external force-of-infection at start of burn-in period 2 years 2 years N/A External force-of-infection (foi) during the simulation Varied between runs, as specified in figure S1; assumed to remain constant during the entire simulated period Varied between runs, as specified in figure S1; assumed to remain constant during the entire simulated period N/A Average mosquito biting rate for adult men (mbr), for a relative biting rate of 1 Varied between runs, as specified in figure S1 Varied between runs, as specified in figure S1 N/A Probability distribution describing variation in the individual exposure index, due to personal factors (fixed through life) given age and sex Varied between runs, as specified in figure S1 Varied between runs, as specified in figure S1 Gamma distribution is assumed as in 2  The total mf density in the population contributes towards the current density of L3 larvae in the human-biting mosquito population, where the distribution of L3 amongst the human-biting mosquito population is completely homogeneous. An empirically derived relationship is used for the uptake of mf by a mosquito, where both Culex and Anopheles uptake curves are implemented depending on setting (see Irvine et al. 1 ). The model dynamics are therefore divided into the individual human dynamics, including age and turnover; worm dynamics inside the host; microfilariae dynamics inside the host and larvae dynamics inside the mosquito.

Model implementation
We performed simulations for LF endemic settings in India and Africa, accounting for different population sizes, possible range of precontrol mf prevalence levels (prevalence range 0 to 60% in the scenarios considered) and regional differences in local vector species and standard treatment regimens (Culex and DA in India, Anopheles and IA in Africa). More specifically, in order to generate the required range of mf prevalences in individuals above 5 years of age, we varied three parameters of the model, the vector to host ratio (V/H), the average population bite risk (k) and the importation rate (Imp), using parameter sets from a range of plausible values based on previously analysed data 1,2,3 . The graphical representation of the values is shown in Figures S6-S8.
For stochastic models it is essential that an importation rate is included, otherwise the equilibrium distribution (steady state) that is used as the starting point of the simulations can potentially converge to the degenerate distribution where no-one is infected. The importation rate does not need to be large, in fact it should not be driving the infection. For this LF study we used a random number drawn from a uniform distribution with minimum 0 and maximum 0.00025 (max 2.5/10000 infections per month). The interventions reduce the prevalence over time, and therefore as years pass, the importation rate decreases in proportion to the reduction 29 in prevalence seen in some pilot simulations. More specifically, we produced 2000 simulations for each scenario with constant importation rate over the years (including the years of MDA). Then, for our main set of simulations, we adjust the importation rate according to how the prevalence changed in our pilot runs after the intervention was applied (for example see Figure  S9).
Finally, compliance between rounds of MDA is modelled based on the paper by Griffin et al. 14 , following the description in Dyson et al. 15 and previous implementation of the model 2 . A summary of all model parameters is available in Table S5. Figure S6 -Vector to Host ratio against aggregation parameter k. The density plot indicates the parameter space areas from the simulations that were in the 0 to 60% mf prevalence range at baseline. Dark blue areas denote more common parameter values.

Figure S7
-Importation rate against aggregation parameter k. The density plot indicates the parameter space areas from the simulations that were in the 0 to 60% mf prevalence range at baseline. Dark blue areas denote more common parameter values.

Figure S8
-Importation rate against Vector to Host ratio. The density plot indicates the parameter space areas from the simulations that were in the 0 to 60% mf prevalence range at baseline. Dark blue areas denote more common parameter values.

Figure S9
-Reduction in mf prevalence over the years seen in the pilot simulations. The first intervention took place in 2018, annual MDA with 65% coverage (IA green and DA red), and was applied for 15 rounds with no interruptions. Proportion of L3 entering host that develop into adult worms 0.00275 [9,10] μ Death rate of adult worms 0.0104 per month [11] δ Production rate of mf per worm 0.2 per month [7] ζ Death rate of mf 0.1 per month [7,12] g Proportion of mosquitoes which pick up infection when biting an infected host 0.37 [13] σ Death rate of mosquitoes 5 per month [8] k Aggregation parameter of individual exposure to mosquitoes Varied Input h(α) Parameter to adjust rate at which individuals of age α are bitten Linear from 0 to 10, with maximum of 1 [9] Imp Importation rate Varied Input