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.
Vector-borne diseases such as malaria are likely to spread with climate change , 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  and vectorial roles . 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 , 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 .
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  to illustrate how temperature may influence disease transmission through vector and parasite development and vector survival.
2. Material and methods
(a) 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.
(b) Data sources for environmental variables and onchocerciasis vectors in forest sites beside the St. Paul river, Liberia
Garms  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 , while vectors in some of its tributaries are S. yahense . 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.
(c) 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 . 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  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).
(d) 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 . 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 , 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 .
(e) 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  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 , 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.
(f) 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 ). 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).
3. Effects of temperature
(a) The impact of temperature on development of Simulium spp.
Cheke  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  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.
(b) 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 .
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' . 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 , larval development and survival times did not depend on how many microfilariae were ingested , 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)  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): 3.1 3.2 3.3where 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.
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  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.  quoted first appearance at 5 days and the MLE predicted 5.8 days. Matsuo et al.  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.  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.  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.
(c) 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.
(d) The impacts of temperature and rainfall on the mortality of the simuliid vector
Le Berre  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  , 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 ). 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).
4. A mathematical model of Simulium damnosum s.l. population dynamics
(a) The model
The following system of ordinary differential equations (ODEs; modified from ), 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). 4.1 4.2 4.3 4.4 4.5where β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 ; μ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.  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 , 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, 4.6 4.7where ɛ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).
(b) Model parametrization: fly fecundity
The fecundity of nulliparous flies is greater than that of parous flies and both vary with fly size . The average sizes of different cytoforms at a given time and place also vary depending on rainfall  via its effect on river discharges . Thus, Cheke & Garms  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  so we used the S. squamosum relations for S. soubrense in Liberia. For this paper, the effect of intraspecific fly size variability is ignored.
(c) Model parametrization: mortality rates
Elsen  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 , but the functional form of that relationship is unknown. Following , 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 , 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: 4.8
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), 4.9where ω = 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 . For S. soubrense in Liberia h = 0.5 based on zoophagy data [13,15] and H = 100 in Gengema . 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).
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.  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.
Edwards & Trenholme  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.
Fly mortality rates were estimated as described in §3(d) and followed a parabolic functional form with temperature (figure 3).
(d) 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, .
(e) 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, .
5. 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  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: 5.1 5.2 5.3 5.4 5.5 5.6where W and M denote the mean number of adult worms per host and of microfilariae per milligram of skin, respectively; L1+2N and L1+2P are the mean numbers of parasite larvae developing in the thoracic muscles of nulliparous and parous flies, respectively (with L1+2P1 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 ); σ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 , 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 .
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 ), 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.
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 , 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 .
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 . 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 . 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.
The views expressed are those of the authors and do not necessarily represent those of DfID, IDRC, the Government of Liberia, UNFCCC or GEF-UNEP.
The part of this study dealing with Ghana was conducted as part of a project entitled ‘Ecohealth approach to the control of onchocerciasis in the Volta Basin of Ghana’ (no. 104270–017), supported by the Climate Change Adaptation in Africa (CCAA) programme, a joint initiative of Canada's International Development Research Centre (IDRC) and the UK Department for International Development (DfID). The Liberian climate change work was part of a study to prepare Liberia's first Climate Change and Vulnerability Assessment National Communication to the United Nations Framework Convention on Climate Change (UNFCCC), supported by the Global Environment Fund of the United Nations Environment Programme (GEF-UNEP). M.G.B., R.A.C., P.H.L.L., M.Y.O.-A. and M.D.W. acknowledge funding from the Wellcome Trust (grant nos. 085133/Z/08/Z and 092677/Z/10/Z), and M.G.B., M.Y.O.-A. and R.A.C. from a Leverhulme Trust–Royal Society Capacity Building Africa Award.
P.H.L.L. is an Imperial College London Junior Research Fellow. M.P. conducted work for this paper as part of her MSc dissertation in epidemiology at Imperial College London. We thank B. A. Boatin, former Director of the WHO Onchocerciasis Control Programme, for supplying some of the data from Ghana.
One contribution of 14 to a theme issue ‘Climate change and vector-borne diseases of humans’.
© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.