The burden of skin disease and eye disease due to onchocerciasis in countries formerly under the African Programme for Onchocerciasis Control mandate for 1990, 2020, and 2030

Background Onchocerciasis (“river blindness”) can cause severe morbidity, including vision loss and various skin manifestations, and is targeted for elimination using ivermectin mass drug administration (MDA). We calculated the number of people with Onchocerca volvulus infection and onchocercal skin and eye disease as well as disability-adjusted life years (DALYs) lost from 1990 through to 2030 in areas formerly covered by the African Programme for Onchocerciasis Control. Methods Per MDA implementation unit, we collated data on the pre-control distribution of microfilariae (mf) prevalence and the history of control. Next, we predicted trends in infection and morbidity over time using the ONCHOSIM simulation model. DALY estimates were calculated using disability weights from the Global Burden of Disease Study. Results In 1990, prior to MDA implementation, the total population at risk was 79.8 million with 26.0 million (32.5%) mf-positive individuals, of whom 17.5 million (21.9%) had some form of onchocercal skin or eye disease (2.5 million DALYs lost). By 2030, the total population was predicted to increase to 236.1 million, while the number of mf-positive cases (about 6.8 million, 2.9%), people with skin or eye morbidity (4.2 million, 1.8%), and DALYs lost (0.7 million) were predicted to decline. Conclusions MDA has had a remarkable impact on the onchocerciasis burden in countries previously under the APOC mandate. In the few countries where we predict continued transmission between now and 2030, intensified MDA could be combined with local vector control efforts, or the introduction of new drugs for mopping up residual cases of infection and morbidity.


Definitions
. List of terminology and definitions used in the manuscript.

Terminology Definition Mf-positive person
Someone who is positive for Onchocerca volvulus microfilariae (mf) in a single skin snip test.

APOC project
A geographical implementation unit for community-directed treatment with ivermectin (CDTi) under the formerly African Programme for Onchocerciasis Control (APOC). Each APOC project has its own organisational structure responsible for implementing the recommended CTDi strategy [1]; a list of all APOC projects is available in Supplement S3.

Onchocercal eye disease (OED)
Functional visual impairment or blindness. Following the WHO criteria, we define visual impairment as visual acuity in the range from worse than 6/18 to 6/60 and equal to or better than 3/60 in the better eye, which covers both moderate and severe visual impairment. According to the same guidelines, blindness was defined as visual acuity of less than 3/60 or a restriction of visual field to less than 10° in the better eye [2].

Onchocercal skin disease (OSD)
We considered six different subtypes of onchocercal skin disease (OSD), namely severe itch, reactive skin disease, palpable nodules, hanging groin, atrophy (<50 years of age), mild and severe depigmentation. We applied a publicly available clinical classification and grading system [3] to define the conditions.

P5-project
A combination of all areas in a country that are considered to be hypoendemic for onchocerciasis (based on REMO surveys showing nodule prevalence levels between 5% and 20%), that were not part of previously defined APOC projects, and that are likely to require treatment or other interventions for the purpose of onchocerciasis elimination [4].

P20-project
A combination of all areas in a country that are considered to be hyperand mesoendemic for onchocerciasis (based on REMO surveys showing nodule prevalence levels of ≥20%), that were not part of previously defined APOC projects, and that are targeted for MDA for the purpose of onchocerciasis elimination [4].

Pre-control infection levels
We

Population size
The population size per APOC project for 1995 was previously collected through a census conducted by community drug distributors for estimating the amount of ivermectin required in mass treatments. For projects where such information was not available, we used the population sizes as previously estimated and published by others [4,8]. We were interested in the population size for 1990 as some APOC projects already initiated treatment between 1990 and 1995. To correct the population size from 1995 to 1990, we used data from the United Nations Population Division [9].
On the basis of these data, we assumed an annual population growth of 2.75%. The pre-control population size for the year 1990 per APOC project was calculated by multiplying (1/1.0275)^5 with the available population size in 1995 (reduction in population size between 1995 and 1990). We estimated the population density for 2020 and 2030 by multiplying the pre-control country-specific population sizes by an assumed annual population growth factor of 1.0275^t, where t is the time in years.

History of control
We used a previously published APOC treatment database which contains information on pre-control endemicity, MDA start year, MDA frequency per year, and treatment coverage per APOC project up to November 2013 [8]. We adopted a recently updated version of the MDA data which includes some additional projects that were previously categorised as not requiring MDA, but do require MDA after all as described elsewhere [4]. For APOC projects where MDA was expected to start between 2013 and 2017, we consulted the ESPEN portal and corrected MDA start year where possible [10]. We also considered ivermectin distribution for lymphatic filariasis (LF) in areas that were endemic for both LF and onchocerciasis. For areas that initiated MDA before 2017 according to ESPEN, we assumed that MDA will continue until at least 2030 at the same coverage and frequency as previously reported [4,8]. The WHO recommended in 2020 to suspend all epidemiological surveys and MDA activities for Neglected Tropical Diseases tackled by preventive chemotherapy and transmission control (PCT) due to severe acute respiratory syndrome coronavirus 2 disease 2019 (COVID-19) pandemic [11].
This has led to severe disruptions of MDA programmes throughout all onchocerciasis-endemic countries. To account for the COVID-19 pandemic in our analysis, we have set the therapeutic coverage of MDA in the year 2020 to zero for all APOC projects; in case bi-annual MDA was provided both rounds were assumed to be missed.
We further assumed that untreated hypoendemic areas suspected of loiasis co-endemicity ("P5 areas", Table A) will remain untreated until at least 2025. This is in accordance with the Mectizan Expert Committee/Technical Consultative Committee (MEC/TCC) guidelines that recommended modified individual treatment rather than MDA in areas hypoendemic for onchocerciasis with known or suspected loiasis [12]. Hypoendemic areas that are non-endemic for loiasis were assumed to start MDA by 2023 with population coverage levels as reported by Kim et al. [8]. The underlying assumption is that L. loa-free hypoendemic areas will start additional mapping of onchocerciasisendemic areas earlier than loiasis-endemic areas and will be able to start MDA more rapidly due to absence of implementation obstacles related to loiasis. reported coverages of >100%. [14] The treatment history in Liberia, which was not included in the CDTi database by Kim et al. [8] and  [16]. ESPEN reports a country-wide therapeutic coverage of 49% in 2015, so we adapted this coverage also to 2014 to account for the EVD outbreak.
We applied the reported MDA epidemiological coverages between 2015 and 2019, as reported by ESPEN. MDA coverage of ivermectin distribution was averaged, as described above. We assumed that all NTD programmes were suspended in 2020 due to the COVID-19 pandemic, including all MDA activities. We further assume that the MDA coverage for 2021 to 2030 remains on average stable similar to the reported coverage of 2019. REMO surveys estimated that the population at risk for onchocerciasis in 1999 was 1.1 million [15], and this number was used to estimate the population at risk in 1995 using the UN population division growth factor [9].
In addition, in the current study we also consider that vector control was systematically implemented on the island of Bioko (Equatorial Guinea) [17,18]

Model background
We used the mathematical model ONCHOSIM to predict trends over time in the prevalence of infection and a wide range of clinical manifestations due to onchocerciasis. ONCHOSIM is a wellestablished mathematical model for simulating transmission and control of onchocerciasis in a dynamic population [26]. A detailed formal description of ONCHOSIM has been published elsewhere [27]. As previously, individual host participation in preventive chemotherapy was assumed to be a mix of random and systematic participation (some people are more inclined to participate than others). Ivermectin is assumed to clear 100% of the mf from the host, and to permanently reduce the capacity of female worms to produce mf by, on average, 34.9% [25]. We further assumed that treatment temporarily stops mf production by adult female worms altogether; mf production then gradually recovers over a period of 11 months on average. As in previous modelling exercises [26,28,29], we assumed no effect of ivermectin on prepatent worms.

Generic disease module
A new model extension for simulating morbidity due to onchocerciasis in ONCHOSIM (WORMSIM version 2.76) was recently developed and quantified, and is described elsewhere [27]. In short, the Together, these three parameters determine the shape of the community-level association between prevalence of infection and a specific condition or continuum of clinical manifestations, as well as age patterns in morbidity at given levels of infection in the community. More information about the biological assumptions and implementation of these disease parameters in the model can be found elsewhere [27].
The pre-control association between the age-patterns in prevalence of each subtype of onchocercal skin disease (OSD) was quantified using data from a large multi-country dataset originating from forest areas [27]. Palpable nodules, severe itch and reactive skin disease (RSD) were considered to be acute, reversible clinical manifestations where regression of symptoms can occur after ivermectin treatment. The speed of reversion was calibrated using longitudinal trends in the prevalence of infection and morbidity after initiation of MDA [30,31]. The pre-control prevalence of onchocercal eye disease (OED) was quantified using various datasets on the community-level prevalence of infection and morbidity from forest and savanna areas [32][33][34][35][36][37][38][39]. Both stages of vision loss (visual impairment + blindness) were considered to be irreversible. More details on parameter quantifications for morbidity can be found elsewhere [27].

Simulations per APOC project
We divided each APOC project into six endemicity categories (Tables A and B) and for each category calculated the mean mf prevalence across pixels [6]. Next, in ONCHOSIM we calibrated the parameter for relative biting rate (rbr) to reproduce the mean pre-control mf prevalence in each project and endemicity category. As in previous ONCHOSIM modelling studies [8,26,29,40], exposure heterogeneity was set at k = 3.5. The rbr parameter was calibrated by simulating pre-control mf prevalences (all ages) for a grid of rbr values with steps of 0.01 (Fig A). Then, for each project and endemicity category, we mapped the mean pixel-level mf prevalence to a value for rbr, using linear interpolation. The corresponding rbr values were used in further simulations for the effect of MDA on infection and morbidity.
S3 Text provides an overview of the history of control and assumptions about future MDA used in the simulations. We assumed that vector control reduced vector abundance by 97.5%. Due to the absence of pre-control pixel-level data from the island of Bioko [6], based on literature we assumed that the overall crude pre-control mf prevalence on the island was 75.2% [41]. To reproduce this very high mean pre-control mf prevalence, we assumed that the island population was distributed over the two highest endemicity categories (with average mf prevalence of 68.5% and 79.1%, respectively) at a ratio of 25:43 such that average prevalence over the whole population was 75.2%.

Fig A.
Association between relative biting rate (rbr, x-axis) and the pre-control mf prevalence (all ages, y-axis) as predicted by ONCHOSIM. The level of exposure heterogeneity was set at k = 3.5.
As ONCHOSIM cannot simulate stable infection with a community mf prevalence of below ~30% (these simulations result in spontaneous fade-out of transmission), we calculated the manifestations of disease prevalence in hypoendemic areas (endemicity level 1) as a ratio of the model-predicted morbidity prevalence in mesoendemic areas (endemicity level 2) ( Table B), as previously done [40].
To assess which ratio in morbidity prevalence between mesoendemic versus hypoendemic areas would correspond to what is reported in the field, we performed a meta-analysis of published data on the prevalence of various subtypes of skin manifestations in hypoendemic versus mesoendemic areas [42,43]. For this purpose, we classified the different subtypes of skin manifestations in three groups corresponding to similar levels of acuteness and reversibility: (1) acute, reversible conditions (i.e. RSD and severe itch); (2) palpable nodules (i.e. slowly reversible); and (3) chronic, irreversible conditions (i.e. hanging groin, atrophy, depigmentation). In the meta-analysis, we estimated the relative difference in prevalence of aforementioned three groups of skin manifestations (i.e. using a log link function), allowing for random effects to capture heterogeneity between the two studies as well as between individual types of skin manifestations within each of the three groups. The results of the meta-analysis are visualised in Fig B. Trends in mf prevalence in hypoendemic areas were derived by scaling trends for mesoendemic areas by the same amount as acute, reversible conditions, as these two metrics follow each other closely over time [27]. in mesoendemic areas (level 2) for these six APOC projects was further disregarded.

Fig B.
Results of the meta-analysis of published data [42,43] on the prevalence of OSD in hypoendemic versus mesoendemic areas. The risk ratio (RR) represents the factor that should be multiplied by the morbidity prevalence in mesoendemic areas for each corresponding clinical manifestation to calculate the morbidity prevalence in hypoendemic areas.

Disease burden calculation
In this section, we describe how the disease burden of onchocerciasis was calculated in terms of disability-adjusted life years (DALYs), which are defined as the sum of Years Lived with Disability (YLDs) and Years of Life Lost (YLLs) [44,45]. DALYs measure the time lost due to the effects of a condition in terms of (1) the time spent disabled by a condition (YLDs), weighted by the severity of the condition; and (2) time lost due to premature mortality by the condition (YLLs). One DALY represents one year of healthy life lost [44,45].

Years Lived with Disability
YLDs were calculated by multiplying the model-predicted number of prevalent cases with symptoms with a weighted severity level of the condition. These so-called disability weights have been developed for a very wide range of conditions as part of the Global Burden of Disease (GBD) studies [46,47]. The disability weights are estimated on the scoring of the severity of a health state by a participant assessed through face-to-face, telephone and online surveys from nine countries across all continents worldwide [44,48]. Although these is some variability in the scoring of disability weights for vision loss and blindness between respondents due to various factors [49], we have used the disability weights for visual impairment and blindness as most recently published by the GBD [44]. As these conditions were designed to be generic and not necessarily specific to any particular disease, we mapped each of the symptoms considered in this study to a condition associated with a disability weight in the GBD project. For OED, we adopted disability weights for functional vision loss directly from GBD (Table C). For OSD, we used a set of disability weights that were previously compiled by the consulting experts on onchocerciasis for the GBD 2010 study (which included coauthors Murdoch, Stolk, and Coffeng). This previously developed scheme considered three general severity levels of OSD (OSD level 1: slight, visible physical deformity that is sometimes sore or itchy, OSD level 2: visible physical deformity that is sore and itchy, OSD level 3: obvious physical deformity that is very painful and itchy), each of which could be associated with itch or not. Where previously, each of aforementioned severity levels encompassed a suit of OSD types, for our current analysis, we mapped each specific subtype of OSD to one of the severity levels. The disability weights and lay descriptions for each of the OSD types are shown in Table C. Person is completely blind, which causes great difficulty in some daily activities, worry and anxiety, and great difficulty going outside the home without assistance.
* RSD was defined as the presence of acute papular onchodermatitis (APOD), chronic papular onchodermatitis (CPOD), and/or lichenified onchodermatitis (LOD). To account for the resulting mix of severities of reactive skin conditions, we calculated an average disability weight, weighted by the prevalence of the three conditions in the multi-country dataset (APOD: 7.44%, CPOD: 12.89%, LOD: 1.23%) [43]. The three disability weights used for APOC, CPOD, and LOD were based on previous discussion with experts (dr. Michele Murdoch, professor Roderick Hay) for estimates by the Global Burden of Disease. APOD was assigned the disability weight for disfigurement: level 1 and CPOD and LOD were assigned the disability weight for disfigurement: level 2.

Concurrence correction of disability weights
To correct YLD estimates for potential concurrence of multiple symptoms in the same individual (which can hypothetically lead to the sum of disability weights for concurring symptoms to exceed 1.0), we apply a multiplicative correction [50] to avoid overestimation of YLDs due to onchocerciasis [51]. Previous work has shown that there was considerable overlap of onchocercal morbidity within individuals from a hyperendemic rainforest area of Cameroon, i.e. predominantly the occurrence of RSD and severe itch [50]. Different subtypes of OSD may lead to similar expressions of disability, e.g.
stigmatisation (such as for depigmentation, hanging groin, atrophy, RSD) or low self-esteem due to skin disfigurement (e.g. RSD, open wounds due to scratching as a result of severe itching) [52]. It has therefore already been suggested by others that any overlap between clinical manifestations due to the same disease should be taken into account in burden estimates [50,[53][54][55]. We choose to add up DALY estimates for OED and OSD without concurrence correction based on the notion that the two groups of symptoms involve completely different mechanisms of imposing a burden on an individual [50]. Any resulting potential overestimation of YLD will be limited because OED is relative rare compared to OSD (i.e. the concurrence correction would only affect a small number of people).
As such, we used a multiplicative approach in the estimation of disability weights to correct for nonrandom overlap of subtypes of OSD in the calculation of disease burden due to onchocerciasis. To do this, we first simulated all possible combinations of subtypes of OSD using ONCHOSIM, after which we obtained the prevalence of condition A and / or B (A|B), A|C, B|C, and A|B|C (here an example of three conditions that potentially concur). We then calculated the number of cases by multiplying the predicted prevalence with the population at risk. We could then calculate the number of cases with condition A and B (A₼B), A₼C, B₼C, and A₼B₼C (see Fig C). For example, condition A₼B was calculated by: Condition A₼B = +condition A(unc) + condition B(unc)1 − condition A|B Where:

Condition A₼B = Persons with both conditions A and condition B
Condition A(unc) = Persons with at least condition A, and potentially also condition B and/or C.
Condition B(unc) = Persons with at least condition B, and potentially also condition A and/or C.

Condition A|B = Persons with condition A and/or B.
In order to calculate the disability weight for each of the combinations of subtypes of OSD possible, we used the disability weights as assigned to each clinical manifestation ( Table C) and applied the multiplicative model [53,54], using the following formula (here example of two concurrent conditions): Where:

DwA₼B =
The disability weight of the combined presence of conditions A and B (A₼B).

DwA =
The disability weight of condition A as provided in Table C.

DwB =
The disability weight of condition B as provided in Table C.
The burden was then calculated by multiplying the corrected disability weight for each possible comorbidity with the number of cases of the respective concurring conditions (here an example of two concurrent conditions):

Burden(A₼B) = DwA₼B * Cases condition A₼B
Where: Burden (A₼B) = the total burden of condition A and condition B, for which it is as yet unknown which proportion of the burden is attributable to condition A and which proportion is attributable to condition B.
As we were interested in the burden of each unique clinical manifestation separately, we calculated the relative proportion of the burden attributable to each condition (here example: DALYA and DALYB). This was done using the following equation (here an example of two concurrent conditions): Relative proportion of YLD @ out of the YLD @₼B = YLD @₼B * DE@ DE@ F DEB and Relative proportion of YLD B out of the YLD @₼B = YLD @₼B * DEB DE@ F DEB The total DALYs for each clinical manifestation corrected for co-morbidity was then calculated by summing for each possible co-morbidity the relative proportion of DALYs due to that subtype of OSD.
Based on a theoretical maximum correction approach in case of maximum non-random concurrence of the six skin manifestations (i.e. severe itch, RSD, depigmentation, atrophy, hanging groin, palpable nodules), we calculated a reduction of 6.7% in total YLDs for OSD as compared to the uncorrected YLDs. To avoid the more than ten-fold increase in number of simulations required to perform the concurrence for each project and endemicity category (i.e. to simulate all symptom combinations A₼B₼ …), and yet be conservative, we simply applied the aforementioned maximum concurrence correction for OSD of 6.7% to all projects and endemicity categories.

Years of Life Lost
The YLL measure captures the years of life lost due to premature mortality compared to the remaining life expectancy based on a counterfactual life table for a person of the same age but without the condition [56]. Of all the clinical conditions considered in this analysis, blindness is the only condition that has been associated with excess mortality [57,58]. We therefore calculated the YLLs for blindness only, stratifying by bioclimate (savanna vs forest areas) and the six endemicity levels.  [59]. This life table was also used to simulate human demography in ONCHOSIM (Fig D). For a life

Sensitivity analysis
We performed univariate sensitivity analyses to assess the impact of biological and programmatic assumptions used in our analysis on the estimated number of cases with skin mf and onchocercal morbidity by 2030. An overview of the assumptions (and alternatives) assessed in these sensitivity analyses is presented in Table E.  between mesoendemic versus hypoendemic areas (Fig E-F Caption: *Alternative ratio of morbidity prevalence in hypoendemic versus mesoendemic areas (see Table E).
The case estimates for OED by 2030 greatly depend on the assumptions applied to the possibility of regression of eye damage once symptoms have developed (Fig G). If we assume a 1% rate of damage regression that allows for healing of eye symptoms with a constant fraction of damage that is healed with every monthly time step, then we predict 358.4 thousand fewer cases with visual impairment or blindness. If we assume higher (60%) excess mortality due to blindness, we predict 26.1 thousand fewer OED cases. As many pre-control hypoendemic areas are in the savanna bioclimate, alternative assumptions on the prevalence ratio among hypoendemic versus mesoendemic areas impact the OED case numbers considerably (191.4 thousand fewer cases).