Potential effects of warmer worms and vectors on onchocerciasis transmission in West Africa.

Development times of eggs, larvae and pupae of vectors of onchocerciasis (Simulium spp.) and of Onchocerca volvulus larvae within the adult females of the vectors decrease with increasing temperature. At and above 25°C, the parasite could reach its infective stage in less than 7 days when vectors could transmit after only two gonotrophic cycles. After incorporating exponential functions for vector development into a novel blackfly population model, it was predicted that fly numbers in Liberia and Ghana would peak at air temperatures of 29°C and 34°C, about 3°C and 7°C above current monthly averages, respectively; parous rates of forest flies (Liberia) would peak at 29°C and of savannah flies (Ghana) at 30°C. Small temperature increases (less than 2°C) might lead to changes in geographical distributions of different vector taxa. When the new model was linked to an existing framework for the population dynamics of onchocerciasis in humans and vectors, transmission rates and worm loads were projected to increase with temperature to at least 33°C. By contrast, analyses of field data on forest flies in Liberia and savannah flies in Ghana, in relation to regional climate change predictions, suggested, on the basis of simple regressions, that 13-41% decreases in fly numbers would be expected between the present and before 2040. Further research is needed to reconcile these conflicting conclusions.

Introduction

Vector-borne diseases such as malaria are likely to spread with climate change [1], but little attention has been paid to how future climatic regimes may affect onchocerciasis, a debilitating disease occurring in sub-Saharan Africa, Central and South America, and the Yemen. Onchocerciasis, or ‘river blindness’, owing to infection with the nematode parasite Onchocerca volvulus and transmitted by blackflies (Simulium spp.), causes visual impairment, blindness, a range of skin lesions and excess mortality [2,3]. It has been estimated that in Africa 37 million people were infected prior to the inception of the Onchocerciasis Control Programme in West Africa (OCP) and the African Programme for Onchocerciasis Control (APOC) [4,5]. In West Africa, the vectors are various cytoforms of the Simulium damnosum complex which differ in their ecologies [6] and vectorial roles [7]. Of the principal vectors in West Africa, S. sanctipauli, S. soubrense and S. yahense are found mostly in forests; S. squamosum principally occurs in highland zones, while S. damnosum s.str. and S. sirbanum are more widespread, but this latter pair are the only common species found in northern savannah zones. Simulium soubrense, S. yahense, S. damnosum s.str. and S. sirbanum occur in Liberia, the source of the forest data presented here; and at least six different cytospecies occur in Ghana [8], the source of the contrasting savannah data analysed in this paper. This diversity means that generalizations about effects of climate change on onchocerciasis transmission require caution unless the particular vector or vectors involved are specified, with similar caveats necessary for different river sizes and bioclimatic zones.

The immature stages (eggs, larvae and pupae) of blackflies (Diptera: Simuliidae) are found in fast flowing, highly oxygenated, water. Assuming an adequate food supply in unpolluted rivers, there is a variety of factors that determine: (i) the rate of development from egg to adult, which is principally governed by temperature; (ii) the densities of populations, which are affected by the development rates, river discharges, adult survival rates, immigration and emigration, and (iii) the geographical range, which depends on habitat type, temperature and river structure. In addition, fly size, fly fecundity and water quality will be influential in determining population densities. Generalizations about S. damnosum complex population dynamics are difficult given the importance of local conditions, particularly river topographies and the vegetation in and around river beds. Rising or falling river heights and discharges may either increase or decrease blackfly breeding opportunities. For instance, excellent breeding sites may disappear when a river floods or be created when previously dry rocks or vegetation become partially submerged leading to the formation of rapids. Furthermore, temperatures vary along rivers, increasing with distances from their sources and between neighbouring rivers [9].

The dependence of onchocerciasis vectors on temperature, rainfall and thus river discharges means that they will be affected by changes in climate. Here, we first examine empirical relationships on fly numbers and environmental variables. Next, we show the effects of temperature on vector development rates and survival, and on development of the parasite within the vector which is also temperature dependent [10,11]. We then developed a novel onchocerciasis vector population model to examine how changes in temperature affect fly numbers. The model was parametrized from data obtained from literature reviews and unpublished results. Much of the data came from four sites, two in Liberia and two in Ghana, for which climate change output from regional climate models (RCMs), downscaled from global climate change models, was linked to hydrological models for the relevant river basins. Finally, we linked the savannah vector model to a framework for the parasite's dynamics in humans and vectors [12] to illustrate how temperature may influence disease transmission through vector and parasite development and vector survival.

Material and methods

Literature review

A literature review identified publications that contained data on larval development of O. volvulus or survival and/or development of Simulium spp. at different temperatures. In addition to reviewing early and grey literature, electronic searches were conducted in 2013 using the databases PubMED and Web of Knowledge, supplemented by electronic searches of the DIALOG library conducted between the early 1980s and October 2000.

Data sources for environmental variables and onchocerciasis vectors in forest sites beside the St. Paul river, Liberia

Garms [13] studied the biology of S. damnosum s.l. in the Bong Range, a forested zone of Liberia, through which the St. Paul river flows. It is now known that the onchocerciasis vector that breeds in the St. Paul river, a member of the S. sanctipauli sub-complex, is S. soubrense [14], while vectors in some of its tributaries are S. yahense [15]. Here we will only consider populations of S. soubrense, which were never subject to larviciding, studied at two sites beside the St. Paul river: Haindi (6°53′45″ N, 10°22′48″ W), where data were collected weekly from October 1968 until December 1969 inclusive, and Gengema (6°54′6″ N, 10°21′44″ W), where data were collected monthly from February 1969 until February 1971 inclusive. Vector biting rates, parous rates and associated precipitation and temperature data were measured in the field or obtained from local sources.

Climate change predictions for the St. Paul river basin

For the Liberian climate projections, monthly and daily means of precipitation, temperature and potential evapotranspiration for the 1961–1990 period were obtained from the FAO New LocClim software and database [16]. Daily values of precipitation, minimum and maximum temperatures covering the period 1961–2040 were obtained from archives of the AMMA-EMSEMBLES ensemble-based runs for West Africa [17,18]. The period 1961–1990 was considered as the baseline period for this study, while the period 2011–2040 was considered as the period for the future scenario (denoted ‘2020s’). The New LocClim data were derived from observed data for Liberia and were used in this study as the observed data for the ‘baseline scenario’ as actual daily observed data for the baseline period were not available. The ensemble climate data used were from two RCMs (HadRM 3P and REMO) and were based on the IPCC SRES A1B scenario experiment. The HadRM 3P model was forced with the boundary conditions of the HadCM3 Global Climate Model (GCM) while REMO was forced with the boundary conditions of ECHAM5 GCM. The HadRM 3P and REMO models were chosen from 10 RCMs used in the AMMA-ENSEMBLES, because they represent the driest and wettest future climatic conditions that can be expected for the St. Paul river basin under the IPCC A1B scenario.

Biases in projections for the ‘2020s’ by the two RCMs were corrected using the ‘delta’ approach, which has been used extensively for corrections in RCM projections in climate change impact studies (e.g. [19,20]). The approach involved (i) computing monthly means of rainfall and temperature for the baseline period and the ‘2020s’ for the RCM data; (ii) determining future changes in precipitation and temperature by contrasting the monthly means for the ‘2020s’ with those of the baseline determined in (i), and (iii) applying the changes determined in (ii) to observed data to obtain the ‘2020s’ climate data for the impact studies.

The baseline and the ‘2020s’ climate data were used to drive Budyko's water balance model of a river basin [21] to predict likely changes in the ‘2020s’ compared with the baseline. The results showed that in the ‘2020s’, daily temperatures in the St. Paul river basin are expected to rise by 1.1–1.3°C, rainfall to either decrease by 3.6% (prediction based on HadRM 3P) or increase by 2% (prediction based on REMO) and the mean total annual river flow in the St. Paul to reduce by 0.7% (REMO) or 25% (HadRM 3P).

Data sources for environmental variables and onchocerciasis vectors in savannah sites of the Black Volta and Pru river basins, Ghana

Savannah regions of Ghana were part of the OCP, when rivers including the Black Volta and Pru were regularly treated with insecticides to kill the larvae of S. damnosum s.l. between 1975 and 2002 [22]. At the study sites used here, Agborle Kame (08°14′04″ N, 2°12′23″ W) on the Black Volta river and Asubende (08°01′01″ N, 00°58′54″ W) on the Pru river, the vectors are almost exclusively the savannah members of the S. damnosum complex, i.e. S. damnosum s.str. or S. sirbanum [8], for which historical OCP data on fly biting rates, parous rates and river discharges prior to larviciding treatments were available. Rainfall and temperature data were taken from the FAO New LocClim software and database [16].

Climate change predictions for the Black Volta and Pru river basins, Ghana

Daily observed data on rainfall and temperature for the period 1961–1990 for the Black Volta and Pru river basins were used as the ‘baseline scenario’ for determining changes in future climate and river flow in this study. The observed data were obtained mainly from the Ghana Meteorological Agency. For the Black Volta basin, climate data covering the part of the basin in Burkina Faso were obtained from the Direction de la Météorologie Nationale, Burkina Faso. The climate change scenarios used in this study were obtained from the ENSEMBLES project [17] for two RCMs (HadRM 3P and REMO) simulated under the IPCC A1B greenhouse gas emission experiment. The HadRM 3P and REMO models were driven by boundary conditions of the HadCM3 GCM and the ECHAM5 GCM, respectively. As for Liberia (§2c), the HadRM 3P and REMO models were chosen because they provided the two extremes of future projections from 10 RCMs over both the Black Volta and Pru basins. Data for both were obtained for the periods 1961–1990 (baseline) and 2011–2040 (‘2020s’). Biases in the RCM data were corrected using the ‘delta’ approach described in §2c.

The impact of changes in rainfall and temperatures on river flow in the Black Volta and Pru basins was assessed by integrating the climate change scenarios in the Soil and Water Assessment Tool (SWAT) hydrological model [23], which had been adapted to the two basins. Analysis of the mean of the climate data from HadRM 3P and REMO revealed that the mean daily temperature and annual total rainfall for the Black Volta basin will increase in the ‘2020s’ by 1.1°C and 3.1%, respectively, relative to the baseline. For the Pru basin, increases in the daily temperature and annual rainfall are projected to be 0.7°C and 2.3%, respectively. The mean annual river flow for the ‘2020s’, based on the mean of the climate data from the two RCMs, is expected to increase above the baseline value by 1.8% for the Pru basin and 0.8% for the Black Volta basin.

Water temperatures

Studies of the effects of temperature on development times of immature stages of S. damnosum s.l. have been based on water temperatures (table 1 shows temperature tolerances of different members of the S. damnosum complex in West Africa [24,25]), but climate change projections provide information on air temperatures so a relationship between the two was needed for conversions. River surface water temperatures vary with ambient temperature, elevation, canopy cover and river width. For instance, temperatures vary from about 24°C to more than 30°C in the St. Paul river but are much more stable, varying only from 22.5°C to 24.5°C in small streams in the Bong Hills (see fig. 2b of [13]). Although different functions will be needed for different rivers, river water temperatures (Tw) are generally linearly related to ambient air temperatures (T), as illustrated by data for the St. Paul river (figure 1a) which can be modelled by the equation: Tw = 0.9844 T − 1.0352 (R2 = 0.74).

Table 1.

Water temperature (°C) ranges of rivers in which different members of the Simulium damnosum complex have been found breeding in West Africa during wet and dry seasons, from [24,25].

taxon	water temperature range (°C), wet season	water temperature range (°C), dry season	water temperature range (°C), overall
S. squamosum	23–26.5	22–29	22–29 (mean 25)
S. yahense	23–25	24–30	23–30 (mean 25)
S. sanctipauli	26.5–28	29–33	26.5–33 (mean 27)
S. soubrense	24.5–28	30–33	24.5–33
S. damnosum s.str.	25–27	27.5–33	25–33 (mean 27)
S. sirbanum	25–26.5	27–33	25–33 (mean 27)

Figure 1.

(a) The relationship between water temperature (Tw, °C) and air temperature (T, °C) in the St. Paul river, Liberia. Data collected at the same times on each of several dates in 1968–1970 and 1989. The fitted equation is Tw = 0.9844Tw − 1.0352 (R2 = 0.7385). (b) The relationship between MBR of S. damnosum s.l. and average monthly ln (discharge in m3 s−1) of the Black Volta river at Agborle Kame, Ghana. Data for August and October to December 1974 and January to October 1975 inclusive. The fitted equation is MBR = 445.22(ln(discharge))2−2558.1ln(discharge) + 3843.7 (R² = 0.8587). (c) The relationship between average monthly ln (discharge in m3 s−1) of the Pru river and rainfall in the previous month (RF, mm); data from June 1957 to August 1967 inclusive. The fitted equation is ln(discharge) = 3E−05RF2 + 0.0117RF − 0.1953 (R2 = 0.5438). (d) The relationship between average monthly ln(discharge in m3 s−1) of the Black Volta river at Bui and rainfall in the previous month (mm); data for March 1951 to November 1975 inclusive. The fitted equation is ln(discharge) = 0.0135RF2 + 0.2384RF + 50.767 (R2 = 0.6675). (Online version in colour.)

Effects of temperature

The impact of temperature on development of Simulium spp.

Cheke [26] collated data for the development of the different immature stages of the S. damnosum s.l. complex, a dataset which was extended by information from Crisp [27] and re-analysed. The temperature data for egg development ranged from 20°C to 33°C, those for larval development from 20°C to 31.5°C, and data for pupal development ranged from 24°C to 31.5°C; however, only five pupal data points were available. When temperature ranges were quoted in publications, weighted means were used in our analyses. All immature stages (eggs, larvae and pupae) showed a decrease in development time with increasing temperature, although data for low temperatures were unavailable. Functions of the form ΔS(Tw) = a exp(−bTw) were fitted to the data, where ΔS(Tw) is the development time of stage s (in days) as a function of water temperature Tw (°C), and s refers to eggs (E), larvae (L) or pupae (P). Parameter values for each stage are given in figure 2a–c.

Figure 2.

Development times of immature stages of S. damnosum s.l. at different temperatures. (a) Eggs, (b) larvae, (c) pupae (Adapted from [26,27]) and (d) the temperature-dependent development function for O. volvulus (data extracted from articles where experiments were conducted in a variety of onchocerciasis vectors Simulium spp., see the electronic supplementary material, S1). Fitted lines are exponential functions, for which the formulae for a–c are given in table 3 (R2 values = 0.473, 0.496 and 0.391 for a, b and c, respectively). The fitted line for D is Duration to L3 = 49.884e−0.08 (where T = mean air temperature (°C); R2 = 0.695).

Table 3.

Parameter definitions and values for blackfly population dynamics model. Note expressions for durations of immature stages account for conversion of air temperature to water temperature.

symbol	description	value	references
E(t, T)	no. eggs at time t and temperature T
L(t, T)	no. larvae at time t and temperature T
P(t, T)	no. pupae at time t and temperature T
V(t, T)	no. vectors at time t and temperature T	V (t, T) = N(t, T) + Ψ (t, T)
N (t, T)	no. nulliparous flies at time t and temperature T
Ψ (t, T)	no. parous flies at time t and temperature T	Ψ (t, T) = V(t, T)/[exp(μVg)−1]
βN(T)	per nulliparous fly rate of oviposition at temperature T
βP(T)ɛN	per parous fly rate of oviposition
ɛN	no. eggs per nulliparous fly	432 eggs for S. damnosum s.str./S. sirbanum; 492 eggs for S. squamosum	[51,52]
ɛP	no. eggs per parous fly	142 eggs for S. damnosum s.str./S. sirbanum; 215 eggs for S. squamosum	[51,52]
ΔE(TW)	duration of egg stage as a function of water temperature TW	11.493 exp(−0.0701TW)	this paper, figure 2a
ΔL(TW)	duration of larval stage as a function of water temperature TW	87.527 exp(−0.0785TW)	this paper, figure 2b
ΔP(TW)	duration of pupal stage as a function of water temperature TW	20.098 exp(−0.0699TW)	this paper, figure 2c
TW(T)	water temperature as a function of air temperature T	TW = 0.9844 T–1.0352	this paper, figure 1a
	per capita background mortality rate of eggs	0.05 d−1	this paper
αE	density-dependent mortality rate of eggs	1.877 × 10−6 d−1 egg−1 (S. soubrense, Liberia)	this paper
		1.519 × 10−5 d−1 egg−1 (S. damnosum s.str./S. sirbanum, Ghana)
μL	per capita mortality rate of larvae	0.3 d−1	[53]
μP	per capita mortality rate of pupae	0.1 d−1	[54]
μV(T)	per capita mortality rate of vectors at temperature T	0.0462 T2 − 2.671 T + 38.99 d−1 (S. soubrense, Liberia)	this paper
		0.0027 T2 − 0.163 T + 2.602 d−1 (S. damnosum s.str./ S. sirbanum, Ghana)
Ψ	proportion of parous flies	1/[exp(μVg−1)]	this paper
g	length of gonotrophic cycle	3.5 days	[44,46]

The impact of temperature on development of Onchocerca volvulus within the simuliid vector

Although unlikely to affect the life history of adult worms in the homoeothermic human host, climate influences the development rate of the parasite within its vector. For instance, in Simulium woodi the worms take 16 days to reach the infective stage at 18°C but only 4 days at 28°C, the highest temperature at which the flies survived; at temperatures between 14°C and 17°C the larvae developed abnormally [28].

Data from 25 articles, describing 49 studies that gave information on both the development time of O. volvulus in a simuliid vector and the temperatures (16.3–30.5°C) at which the studies were carried out, were analysed. Mean temperatures were used when given, otherwise the midpoint of the provided range was taken as the mean. Times of first appearance of third stage (infective) larvae of O. volvulus (L3) were used whenever possible, but some studies used vague terminology such as ‘the development to abnormal sausage (pre-L3) stages took about 30 days' [28]. If a range was given, the shortest quoted time was used. The times for the first L3 appearances were chosen, because these were the most frequently quoted data and the most biologically meaningful. It also seemed more suitable for comparing data from such a wide variety of sources because selecting the midpoint of the range requires strong assumptions, particularly if development times are heterogeneous. It was assumed that although heavily infected flies may die quicker than flies with few larvae [29], larval development and survival times did not depend on how many microfilariae were ingested [30], so the microfilarial intake was not incorporated in the model. When second stage (L2) or L3 larvae were seen within the first 3 days, these larvae were ignored as being probably from a previous infection (as the fly-feeding experiments were rarely conducted with newly emerged flies, but mostly with wild caught host-seeking flies, which might have previously fed); microfilariae were grouped with the L1 stage when reported in the thorax. When the flies were checked every 24 h and the dead ones dissected to check for the presence of different larval stages [31-42], the developmental rates were estimated by maximum likelihood (using a multinomial distribution) [43] by fitting the following equations to data on larvae in each stage, L1(t), L2(t), L3(t), as a proportion of the total number of parasite larvae observed on each dissection day (t):where l1, l2 and l3 are the proportions of larvae in each stage, and υ1 and υ2 are the rates of progression from L1 to L2, and from L2 to L3, respectively. A model without a delay between the first and the second stage larvae was tried initially, but including a delay gave a significantly better fit (according to the likelihood ratio test). The development rates were then converted into durations of development from L1 to L2, and from L2 to L3 by taking the reciprocal of the respective rates (table 2), and these values summed to obtain the overall duration from L1 to L3.

Table 2.

Parameter values for rates of progression between Onchocerca volvulus first and second stage larvae (υ1), and between second and third (infective) stage larvae (υ2) in simuliid vectors. EIP: extrinsic incubation period (development time between microfilaria and infective larva). Parameters estimated by maximum likelihood.

species and locality	temperature (°C) (range)	delay as L1 (days)	υ1 (d–1)	υ2 (d–1)	EIP (days)	references
Simulium guianense s.l. S. Venezuela	20.0 (16.0–24.0)	5.56	1.413	0.331	9.28	[33]
Simulium damnosum s.l. Cameroon	21.5 (19.5–23.5)	4.37	0.552	0.569	7.94	[37,38]
Simulium ochraceum s.l. Guatemala	25.0	4.54	0.771	0.525	7.75	[31]
Simulium guianense s.l. S. Venezuela	25.0 (22.0–28.0)	3.32	0.455	0.610	7.16	[36]
Simulium damnosum s.l. Nigeria	27.0 (26.0–28.0)	4.03	—	—	5.09	[32]
Simulium exiguum s.l. Ecuador	27.5 (25.0–30.0)	2.53	0.449	0.518	6.69	[35]
Simulium oyapockense s.l./ S. incrustatum	28.5 (27.0–30.0)	3.77	0.742	1.565	5.76	[42]

Generally, the model of equations (3.2) and (3.3) gave good fits except when the transition from one larval stage to the next was highly synchronous. Eichner [37] quoted first appearance of L3 at 6 days post infection and all development complete by 9 days; the maximum likelihood estimation (MLE) gave an estimate of 7.9 days. Grillet et al. [42] quoted first appearance at 5 days and the MLE predicted 5.8 days. Matsuo et al. [31] did not check the flies on day 4 post infection, so it was assumed that the numbers and type of larvae observed on that day were the same as on day 3; the MLE fit estimated first appearance at 7.7 days and the first observation was at 8 days. Takaoka et al. [33] quoted seeing the first L3 on day 9 and the MLE fit predicted 9.3 days, but there were only 34 flies in their study. Basáñez et al. [36] first observed L3 in S. guianense at day 6, for which the MLE fit estimated 7.2 days. Not all published studies contained data suitable for estimation of υ1 and υ2 separately.

The development rates of O. volvulus for all simuliid species (collated in electronic supplementary material, S1) were used to obtain a temperature-dependent larval development function of the form ΔL1+2(T) = a exp(−bT), where ΔL1+2(T) is the duration in days of parasite larval development from L1 to L3 within simuliid vectors as a function of temperature (figure 2d). At temperatures above approximately 25°C, development can take less than 7 days and so, assuming a gonotrophic cycle of 3.5 days and that an infection was picked up at the first bite, a fly could transmit infective larvae as soon as after its second gonotrophic cycle.

The impacts of temperature and rainfall on blackfly populations and projections based on regressions of the mortality of the simuliid vector

Monthly data on mean numbers of flies caught per person, mean percentages of flies classified as parous, maximum, minimum and mean temperatures (°C), rainfall (mm) and river height (m, as measured at gauges in the St. Paul river) were available for Haindi and Gengema (where the vector is S. soubrense) and similar data but with river discharges (m3 s−1) instead of water levels were obtained for Agborle Kame and Asubende (where the vectors are the savannah species S. damnosum s.str. and S. sirbanum). To convert the Liberian water level data into approximate discharges, a relation linking the two for the Amou river in Togo was used, for which: ln(discharge) = 3.3 river height (m) − 1.29 (P < 0.0001; OCP data).

Figure 1b shows that fly numbers at Agborle Kame are dependent on river discharges, with higher numbers at low or high discharges than at intermediate levels. However, temperature is also an important variable as analysis of variance of the Liberian and Ghanaian data combined revealed that effects on monthly biting rates (MBR) of both average temperature (Tav) and ln(discharge) were significant (F = 4.99, P < 0.03 and F = 5.66, P < 0.02, respectively, d.f. = 1/75). Since there was also a significant effect according to forest or savannah vectors (F = 91.46, P ≪ 0.0001), the following equations were derived: for the Liberian sites, ln(MBR) = 20.54–0.438 Tav − 0.0466 ln(discharge) (P < 0.002); and for the Ghanaian sites, ln(MBR) = 12.32–0.197 Tav + 0.152 ln(discharge) (P < 0.006). Analysis of data on the percentages of parous flies revealed a significant effect of vector species (F = 204.96, d.f. 1/75, P ≪ 0.0001), minimum temperature (Tmin, F = 13.06, P = 0.0005) and ln(discharge) (F = 5.90, P < 0.02), yielding the following equations: for the Liberian sites, asin() = 0.03995 − , where ψ is the proportion of parous flies, and for the Ghanaian sites, asin() = .

Using the average monthly temperature (26.5°C) and the average estimated monthly discharge (214.25 m3 s−1) for the period 1973–1975 for Walker Bridge (7°33′ N, 9°53′ W), upstream of Haindi (River Discharge Database at http://www.sage.wisc.edu/riverdata/scripts/station_table.php?qual=32&filenum=2221), the above relation for Liberia predicts an MBR of 5901 flies per person. Assuming a forecast of a 1.2°C increase in temperature and a 25% reduction in river discharges, then only 3536 biting flies per person per month would be expected. Similarly for Ghana, using the mean monthly temperature (26.8°C) and the mean estimated discharge (128.3 m3 s−1) for Agborle Kame in 1975, the above relationship predicts an MBR of 2412 flies per person. Assuming a forecast of a 1.1°C increase in temperature and a 0.8% increase in river discharge, then only 1849 flies per person per month would be expected. A similar exercise for Asubende, using the mean temperature for 1979 (27.9°C) and mean monthly discharges for 1976 (494 m3 s−1) provides an estimate of 2206 bites per person per month. Assuming a temperature increase of 0.7°C and a 1.8% increase in river discharge this would reduce to 1865 bites per person per month. Such predictions could be improved by direct links to rainfall, commonly considered in the output of climate change models, since in both the Pru and Black Volta rivers discharges are related to their basins' composite rainfall in the previous month (figure 1c,d). However, the above extrapolations conflict with results from the dynamic model described in §4(a) and by calculations using matrix population models (RA Cheke 2012, unpublished data), possibly because they do not take into account the feedback loops and regulatory processes incorporated in the population dynamics model (i.e. essentially nonlinear processes); instead, they are based on linearized statistical relationships that do not account for the underlying biology. However, as will be seen in §4(a), the model includes temperature dependency of development and survival rates but does not yet include the effect of rainfall, warranting further work.

The impacts of temperature and rainfall on the mortality of the simuliid vector

Le Berre [44] pointed out that there are three main types of relationships between rivers and blackfly numbers: (i) those in which the flies' numbers are synchronous with river levels; (ii) those in which they vary inversely with river levels and (iii) those in which they are bimodal, such as at Agborle Kame, with peaks related to the availability of larval supports, such as vegetation, being optimal at low water levels and then again at high levels. Given this complexity it is difficult to generalize about effects of rainfall on fly survival or, indeed, on determinants of carrying capacities, so for this paper henceforth we restrict treatment of environmental influences on fly mortalities to temperatures. Furthermore, although there is some agreement that temperatures will increase, climate change projections for precipitation changes in West Africa vary according to the model used, although the likelihood of an increase in extreme events is generally accepted.

The proportion of parous flies (ψ) was used to calculate the probability of daily survival, p, by applying the parous rate formula [45] , where g is the average duration between two consecutive blood meals (i.e. assuming gonotrophic concordance, the mean duration of the gonotrophic cycle), and the daily mortality rate, μV, by taking –ln(p) under the assumption of exponential distribution of survival times and no difference in survival between nulliparous and parous flies. These mortality rates were calculated for the Haindi and Gengema data combined, as well as for the Agborle Kame and Asubende data combined, assuming g = 3.5 days [44,46], and were fitted by least-squares estimation to the monthly average temperature (Tmav) with a polynomial of the form aTmav2 + bTmav + c. Since the parous rates of S. soubrense in Liberia were very low (averaging 15%), the resulting mean life expectancy of the flies was approximately 2 days, possibly highlighting deficiencies with this method. By contrast, the mean parous rate of S. damnosum s.str./S. sirbanum in Ghana was 64%, and the mean life expectancy was 10 days (of the same order of magnitude as that estimated with laboratory data [29]). The parameter values for the flies in Liberia were a = 0.046; b = −2.671 and c = 38.99 (figure 3a); the values for the flies in Ghana were a = 0.003; b = −0.163 and c = 2.602 (figure 3b).

Figure 3.

Mortality rates (μV) of (a) S. soubrense at Gengema and Haindi, Liberia, in relation to average monthly temperature (°C) and (b) for S. damnosum s.str./S. sirbanum at Agborle Kame and Asubende. For means of calculation of the mortality rates from parous rates see text. The equation for the fitted line for (a) is μV = 0.0462Tav2 – 2.671Tav + 38.988 (R2 = 0.426) and for (b) is μV = 0.0027 Tav2 − 0.163 Tav + 2.602 (R2 = 0.031). Data from [13]. (Online version in colour.)

A mathematical model of Simulium damnosum s.l. population dynamics

The model

The following system of ordinary differential equations (ODEs; modified from [47]), describes the rates of change with respect to time (and dependent on temperature) of simuliid eggs E(t, T), larvae L(t, T), pupae P(t, T) and adult females distinguished between nulliparous N(t, T) and parous Ψ(t, T) flies. The total vector population, V(t, T) = N(t, T) + Ψ(t, T).where βN(T) and βP(T) are the per nulliparous and per parous fly daily rates of oviposition (of viable eggs), respectively, which depend on temperature as described in §3(a); 1/ΔE(T) is the rate of progression from eggs to larvae (the reciprocal of the duration in the egg stage, dependent on temperature with a negative exponential relationship as described above); μE[E(T, t)] is the per capita mortality rate of eggs (dependent on egg density according to Kyorku & Raybould [48]; μL and μP are the per capita mortality rates of larvae and pupae, respectively; 1/ΔL(T) is the rate of progression from larvae to pupae; 1/Δp(T) is the rate of eclosion from pupae of adult (nulliparous) flies, half of which are females; μV(T) is the per capita mortality rate of flies (which depends on temperature as described in §3d), and g is the mean duration of the gonotrophic cycle (as in §3d), which according to Takaoka et al. [11] decreases (by up to 2 days) with increasing temperature but does not vary substantially between 18°C and 28°C.

In blackflies, the adult host-seeking female population V(t,T) is often classified as nulliparous (those emerging from pupae and coming to bite for the first time), and parous (those which have fed before, laid eggs and are coming to bite for the second or third time, etc.), but unlike mosquitoes, it is not possible to distinguish between 1-parous, 2-parous, 3-parous, etc. The value of g is usually considered as being between 3 and 4 days [44,46], but may be as low as 2 or 2.5 days [11,49,50]. Following [47], the daily oviposition rates, βN(T) and βP(T), are estimated as the expected number of eggs oviposited in a fly's lifetime divided by the life expectancy. For nulliparous flies, the latter is the reciprocal of the rate at which they progress to the parous compartment (i.e. the duration of the gonotrophic cycle), and for parous flies it is the reciprocal of the parous mortality rate under the assumption of an exponential distribution of survival times. Therefore,where ɛN and ɛP are the average number of (viable) eggs per oviposition laid by nulliparous and parous flies, respectively.

The expression exp(−μV(T)g) + exp(−2μV(T)g) + exp(−3μV(T)g)+⋯ is the proportion of parous flies surviving consecutive gonotrophic cycles with mortality rate μV(T). Nulliparous and parous flies have different fecundity rates (see §4b).

Values used for the parameters are explained in §4(b,c) and given in table 3. The model and the framework described below (§5) were coded in Berkeley Madonna v. 8.0.1 [55] using Runge–Kutta 4.

Model parametrization: fly fecundity

The fecundity of nulliparous flies is greater than that of parous flies and both vary with fly size [51]. The average sizes of different cytoforms at a given time and place also vary depending on rainfall [56] via its effect on river discharges [57]. Thus, Cheke & Garms [52] derived the following expressions for the fecundity (eggs laid per oviposition) for nulliparous S. damnosum s.str./S. sirbanum: no. of oocytes (O) = 1370.79 thorax length (λ) in mm − 939.06, and for parous S. damnosum s.str./S. sirbanum: O(λ) = 1063.95λ − 922.06; for nulliparous S. squamosum: O(λ) = 1230.65λ − 738.69, and for parous S. squamosum: O(λ) = 1081.33λ − 866.77. So, a nulliparous S. squamosum with a thorax length of 1.0 mm would be expected to lay ɛN = 492 eggs and a similar-sized S. damnosum s.str./S. sirbanum nulliparous fly would lay 432 eggs, with parous flies of similar sizes laying ɛP = 215 and 142 eggs, respectively. Data on fecundities of members of the S. sanctipauli sub-complex are lacking but their adults are known to be larger than S. damnosum s.str./S. sirbanum when they are sympatric [56] so we used the S. squamosum relations for S. soubrense in Liberia. For this paper, the effect of intraspecific fly size variability is ignored.

Model parametrization: mortality rates

Eggs

Elsen [58] reported that only 2.7% of eggs reached the pupal stage, which for a development period from egg to pupa of roughly 14 days, yields a daily mortality rate of 0.77; in the absence of reliable field data for egg survival, it is therefore assumed that the per capita mortality rate of eggs is μE = 0.8 d−1. However, given that this was measured in the field, the value corresponds to overall mortality and not to background mortality. In addition, it has been reported that egg mortality increases with egg mass density [48], but the functional form of that relationship is unknown. Following [59], a linear function was chosen, so that with the background rate of egg mortality and αE = the rate of excess mortality per additional egg in the egg mass. Following [47], and to stabilize the population, the mortality rate of eggs was parametrized in terms of a carrying capacity (K) of adult vectors. Therefore, the differential equation for the simuliid eggs can be re-written as follows:

So that . The value of was taken as 0.05 d−1.

To derive an expression for K, the equations for the blackfly population dynamics were set to zero and equilibrium expressions for each stage were obtained (omitting temperature dependence for simplicity; see electronic supplementary material, S2, where the formula for the basic reproduction number or ratio, R0, of the blackfly populations, , based on equations (4.2)–(4.5) and (4.8) is also given; if , the unique positive equilibrium of the model is globally asymptotically stable (electronic supplementary material, S3),where ω = 2ΔE(1 + μLΔL)(1 + μPΔP) (1 + μVg) μV.

In order to obtain values for K from equation (4.9), it was necessary to estimate the equilibrium vector population, V*, which was calculated for the Liberia and Ghana data on pre-control daily vector biting rates (DBR*), human population size (H), length of gonotrophic cycle (g) and proportion of blood meals taken on humans (h), following [12]. For S. soubrense in Liberia h = 0.5 based on zoophagy data [13,15] and H = 100 in Gengema [60]. For S. damnosum s.str./S. sirbanum in Ghana, h = 0.74 (authors' unpublished data on blood-meal analyses providing proportions of flies biting cattle, sheep, goats, pigs, dogs and unknown non-human hosts) and H = 462 (the average of the human population sizes in Asubende and Agborle Kame assessed by Ghana Health Services in 1988).

Larvae

Although no data on larval mortality were available, the per capita larval mortality was assumed to be μL = 0.3 d−1 since Davies et al. [53] found that a survival rate of 0.7387 per day was needed to stabilize their model that used field estimates for most of its other parameters.

Pupae

Edwards & Trenholme [54] reported a proportion of pupal survival of 0.75 over a 3-day experiment, equivalent to 0.91 per day, so a per capita pupal mortality rate of μP = 0.1 d−1 was assumed.

Adults

Fly mortality rates were estimated as described in §3(d) and followed a parabolic functional form with temperature (figure 3).

Model simulations of Simulium soubrense population dynamics in Liberia

The model was run until equilibrium values were obtained for the state variables for varying temperatures using values for K derived from mean biting rates at Haindi and Gengema for the air temperature range of 25–29°C, and calculating pre-imaginal development rates according to water temperatures derived from the regression between air and water temperatures. The model predicted total numbers of biting flies, nulliparous flies and parous flies as depicted in figure 4a. This parametrization of the model only gives positive results between 26°C and 32°C and indicates that S. soubrense populations would peak at 29°C before declining. Over this temperature range, the model predicts parous rates to vary as shown in figure 4b. With T = 30, .

Figure 4.

Output of blackfly population models with parameters for Liberia (forest) (a,b) and Ghana (savannah) (c,d). (a,c) Equilibrium densities of total numbers of flies (diamonds), nulliparous flies (squares) and parous flies (circles) at different temperatures (°C). (b,d) Proportions of parous flies at different temperatures (°C). (Online version in colour.)

Model simulations of Simulium damnosum s.str./Simulium sirbanum population dynamics in Ghana

At equilibrium for a range of temperatures, using values for K derived from average biting rates at Agborle Kame and Asubende for the air temperature range of 25–31°C, and calculating pre-imaginal development rates according to water temperatures as described in §3(a), with parameter values as in table 3 and adult mortalities as described for S. damnosum s.str./S. sirbanum in §3d, the model predicts total numbers of biting flies, nulliparous flies and parous flies as depicted in figure 4c. This parametrization of the model gives positive results above 18°C and shows that populations will peak at 34°C before declining. No simulations were run at temperatures of more than 40°C. The model's parous rate predictions are shown in figure 4d. With T = 30, .

A mathematical model of the transmission dynamics of Onchocerca volvulus

In the light of the above, it is possible to modify the model presented in [12] to account for the population dynamics of vectors and the effects of temperature on simuliids and on O. volvulus within the simuliids. The equations for the modified model are as follows:where W and M denote the mean number of adult worms per host and of microfilariae per milligram of skin, respectively; L1+2 and L1+2 are the mean numbers of parasite larvae developing in the thoracic muscles of nulliparous and parous flies, respectively (with L1+2 being larvae in flies in their first parous cycle), is the mean number of infective larvae in parous flies; H is the human population density; h is the proportion of blood meals taken on humans; is the probability of parasite establishment within the human host, dependent on the intensity of exposure to infective larvae, with equation , omitting the time and temperature dependencies to facilitate notation; δV(M) is the probability of parasite establishment within the vector host, dependent on the intensity of microfilarial infection in the skin, with equation , again omitting time and temperature dependencies, with and representing the maximum probability of parasite establishment within humans and vectors that applies when the intensity of exposure to infective larvae or the intensity of microfilarial infection tend to zero; is the probability of parasite establishment within humans when the intensity of exposure to infective larvae is very large (with , the latter can be set equal to zero), and cH and cV represent the magnitude of (negative) density dependence operating upon parasite establishment within humans and vectors, respectively (the latter can be set equal to zero for forest flies [43]); σW, σM and are the per capita mortality rates of adult worms, microfilariae and infective larvae (ignoring the mortalities of L1 and L2 larvae as it was not possible to estimate these from the data, and thus it is assumed that once in the thorax, parasite larvae will progress to the L3 stage); μH is the mortality rate of humans; μV is the mortality rate of vectors, which in addition to being temperature dependent, is also dependent on the microfilarial load ingested by the flies [29], with functional form , where is the background rate of mortality dependent on temperature as described in §3(d), and αV is the per ingested microfilaria excess mortality rate of the vectors (assuming that this rate is independent of temperature); finally, aH is the proportion of infective larvae shed per bite. Parameter values used are the averages given in table 2 of Basáñez & Boussinesq [12].

Results of running the model with S. damnosum s.str./S. sirbanum (Ghana data) parameters for 100 years to ensure endemic equilibrium at air temperatures from 19°C to 33°C are shown in figure 5. The annual transmission potential (ATP, i.e. the number of infective larvae potentially received annually per person [61]), numbers of female worms per human host and numbers of microfilariae per milligram in the skin of human hosts all gradually rise with temperature, but the rate of increase begins to decelerate above 30°C.

Figure 5.

Output of combination of blackfly population model with parameters for Ghana (savannah) and model for onchocerciasis in the human host showing variation of equilibrium values of numbers of female worms per host (triangles), number of microfilariae per milligram of skin (circles) and ATP of the flies (diamonds) with temperature. (Online version in colour.)

Discussion

The extrapolations based on regressions of field data on blackfly numbers with temperatures and river discharges suggested that fly populations are likely to decrease with increasing temperatures and decreasing, or slightly increasing, river discharges, with the magnitudes of the assumed environmental changes based on regional climate change model scenarios. By contrast, the output of the ODE model, albeit lacking explicit inputs regarding precipitation or river discharges, predicted that fly numbers would have a humped distribution peaking at average air temperatures of 29°C in Liberia and 34°C in Ghana. Given that recent average monthly air temperatures in the study sites were 26.5°C and 26.8°C in Liberia and Ghana, respectively, and are expected to increase by up to 1.1°C or 1.3°C on the basis of conservative climate change scenarios, both countries are likely to see substantial increases in numbers of onchocerciasis vectors in the next few decades on the basis of the ODE model, because of accelerating rates of blackfly development with increasing temperature. However, although there is more uncertainty regarding future changes in precipitation, climate change models for Liberia suggest that the discharge of the St. Paul river will decrease by 0.7–25%, but in Ghana increases of 0.8% in the Black Volta river and 1.8% in the Pru river are expected. Depending on how such changes affect fly numbers it is possible, if river topographies allow lower discharges to lead to fewer flies, that in forested areas of West Africa such as Liberia onchocerciasis transmission may decrease (as implied by the relations based on empirical data) but in savannah zones such as northern Ghana it might increase. The latter conclusion is supported by the results obtained when the vector model was linked to the dynamics of the parasite in humans and vectors which indicated increases in worm burden with increasing temperature in Ghana up to temperatures of 33°C, because of accelerating development of parasite larvae within vectors. Reconciling the contrasting results between those from the field data statistical relationships and the ODE dynamics model is imperative and highlights the deficiencies of statistical (linearized) relationships as opposed to nonlinear dynamics on the one hand, and the need to refine the dynamic model by including dependencies of the vector carrying capacity with rainfall and river levels on the other hand. Also, the modelling of vector mortality was relatively simplistic, assuming an exponential distribution of survival times and a constant (and equal) mortality rate for nulliparous and parous flies. Further work is necessary to better understand the dependency of vector survival and increasing temperature, and laboratory data on neotropical vectors (not shown) suggest that vector mortality may increase with temperature at rates higher than those derived from the observed proportions of parous flies described here. Also, the statistical relationships were derived from existing temperature and river discharge ranges, whereas the dynamic projections predict peak fly numbers to occur above current maximal temperatures, so it is possible that the shape of the empirical functions will change as temperature and rainfall patterns become more extreme.

Although there have been previous studies modelling blackfly population dynamics most have been based on difference equations [53,59,62] and, so far as we are aware, the differential equation model presented here is the first of its kind. While it was designed to be as realistic as possible on the basis of current knowledge, it could be improved by incorporating the effects of rainfall, perhaps by linking precipitation to river discharges, to describe the carrying capacity of eggs explicitly, as opposed to the current construct derived from the mathematical characteristics of the functions used and based on crude estimates of host populations available to the flies, which are probably very low as they are based on sizes of villages where the flies were caught. In addition, more information is needed on the effects of environmental variables on fly mortalities. Indeed the differences in the dynamics of the Liberian and Ghanaian flies may be accounted for by the different mortality/temperature relations used. It is also possible that the gonotrophic cycle length of S. damnosum s.l. is either less than the 3.5 days assumed or that it is temperature dependent, as is the case with other species [11], which will be analysed in future modelling work. Also, the relations for development times of the immature stages of the flies are based on few data and only for water temperatures ranging from 20°C to 33°C, so caution is needed for any forecasts of fly numbers above 33°C. Furthermore as with most, but not all, population dynamics models, immigration and emigration are ignored, even though it is known that savannah members of the S. damnosum complex may migrate up to 300 km [63].

Irrespective of how S. damnosum s.l. population densities will alter with climate change, it is likely that increasing temperatures will lead to changes in the geographical distribution of some species. For instance, S. squamosum and S. yahense are adapted to colder water temperatures than S. damnosum s.str. and S. sirbanum (table 1) so the latter may replace the former at some sites. Such replacements of forest species by savannah species have already occurred in response to habitat changes and are likely to lead to deteriorating epidemiological outcomes [64,65]. Other ecological changes may also have already occurred, perhaps in response to climatic changes during the past 40 or more years. However, even if we had had access to a complete 1974–2001 dataset for the OCP area, any trends would have been difficult to discern there, as they are likely to have been masked by OCP's vector control activities.

The research presented here is timely, as in 2012 the Disease Reference Group for Helminth Infections (DRG4) of the UNICEF/UNDP/World Bank/WHO Special Programme for Research and Training in Tropical Diseases (TDR) identified a need to develop models to investigate the effects of climate change on helminthiases and their control. They recommended conducting literature reviews, experimental/observational studies and parameter estimation to calibrate models on the interaction between the biology of the infections and climate-driven environmental variables [66]. Among the helminth infections transmitted by haematophagous vectors was onchocerciasis, a neglected tropical disease (NTD) earmarked for elimination in the American continent by 2015 and in selected African countries by 2020 according to the World Health Organization (WHO) roadmap for accelerating progress to overcome the impact of NTDs [67]. As our results are inconclusive and conflicting, it is clear that further research is needed before the level of uncertainty surrounding how climate change will affect onchocerciasis transmission, and its control can be reduced. Such further work could encompass seasonality, fly migrations including those of savannah species into forest areas and vice versa, modelling of transmission in the forest (not covered here through lack of space), spatio-temporal changes in human and non-human blood host populations, effects of chemotherapeutic treatments, spatial variation in parous rates and climate-related variations in parameters such as gonotrophic cycle lengths.