Mathematical modeling of climate change and malaria transmission dynamics: a historical review

Abstract

Malaria, one of the greatest historical killers of mankind, continues to claim around half a million lives annually, with almost all deaths occurring in children under the age of five living in tropical Africa. The range of this disease is limited by climate to the warmer regions of the globe, and so anthropogenic global warming (and climate change more broadly) now threatens to alter the geographic area for potential malaria transmission, as both the Plasmodium malaria parasite and Anopheles mosquito vector have highly temperature-dependent lifecycles, while the aquatic immature Anopheles habitats are also strongly dependent upon rainfall and local hydrodynamics. A wide variety of process-based (or mechanistic) mathematical models have thus been proposed for the complex, highly nonlinear weather-driven Anopheles lifecycle and malaria transmission dynamics, but have reached somewhat disparate conclusions as to optimum temperatures for transmission, and the possible effect of increasing temperatures upon (potential) malaria distribution, with some projecting a large increase in the area at risk for malaria, but others predicting primarily a shift in the disease’s geographic range. More generally, both global and local environmental changes drove the initial emergence of P. falciparum as a major human pathogen in tropical Africa some 10,000 years ago, and the disease has a long and deep history through the present. It is the goal of this paper to review major aspects of malaria biology, methods for formalizing these into mathematical forms, uncertainties and controversies in proper modeling methodology, and to provide a timeline of some major modeling efforts from the classical works of Sir Ronald Ross and George Macdonald through recent climate-focused modeling studies. Finally, we attempt to place such mathematical work within a broader historical context for the “million-murdering Death” of malaria.

1 Introduction

Malaria, a potentially deadly disease caused by protozoan parasites known as Plasmodium that infect and replicate within human blood cells, is spread between humans via the bite of the infected female adult Anopheles mosquito, and is one of the greatest infectious maladies to beset mankind. There are five (previously four) Plasmodium species that commonly infect humans, namely P. falciparum, P. vivax, P. ovale, P. malariae, and, very recently, P. knowlesi (Antinori et al. 2012). Of these, P. vivax and P. falciparum are preeminent by far, responsible for nearly all malaria deaths in 2015 (estimated at 438,000 by the World Health Organization (WHO) (WHO 2015), although another major estimate is appreciably higher, at 631,000 deaths (Gething et al. 2016), and the confidence intervals for both estimates are broad). Over 90% of all malarial mortality is attributable to P. falciparum in sub-Saharan Africa, where children under the age of five are chiefly burdened (WHO 2015), and Fig. 1 demonstrates the concentration of malaria risk in this region. Consequently, the focus in this paper is almost exclusively on P. falciparum malaria in Africa.

Fig. 1

Global populations at risk of malaria, in 2013. Tropical Africa is at highest risk, with many countries having 100% of their populations at risk; mortality is also strongly concentrated in this region. Map generated by the World Health Organization’s Malaria Mapper (http://www.worldmalariareport.org/node/68), based on the World Malaria Report, 2015 (WHO 2015)

The emergence of P. falciparum as a major human disease, likely dating back to the acquisition of P. falciparum from a gorilla in Africa some 10,000 years ago (Loy et al. 2017; Carter and Mendis 2002), was directly linked to environmental changes, namely, the end of the last ice age leading to an era of global warming and the subsequent birth of human agricultural civilization, which, via land-use changes and the concentration of human settlement, allowed malaria and its mosquito vectors to thrive (Carter and Mendis 2002; Webb 2014; Packard 2007). The parasites and vectors, having temperature- and rainfall-dependent lifecycles, are restrained by climate to the globe’s warmer latitude and altitude ranges (Patz et al. 1996). Thus, in the modern era, anthropogenic global warming, driven principally by fossil fuel combustion, but secondarily by global land-use changes (IPCC 2013) (primarily deforestation (IPCC 2013), which is in turn driven mainly by agriculture (Rudel et al. 2009; McKinley et al. 2011)), threatens to expand the potential range, and possibly the overall burden as well, of malarial disease. Aside from global restraints, malaria incidence follows altitude in multiple countries, such as Zimbabwe and Kenya (Patz et al. 1996), and the recent expansion of disease into some upland areas, notably the highlands of western Kenya, may be at least partly attributable to warmer temperatures (Pascual et al. 2006; Pascual and Bouma 2009).

However, the ultimate effect of climate change on malaria is far from certain, as a wide milieu of social, biotic, and abiotic factors influence the disease in non-linear ways, and the global burden of malaria contracted enormously over the twentieth century in the face of modest warming (Carter and Mendis 2002; Gething et al. 2010) (although this pattern generalizes poorly to malaria in Africa (Carter and Mendis 2002)). Over the last few decades, a number of mathematical models, typically statistical (using data and statistical approaches to correlate some climate variables with malaria incidence) or mechanistic (accounting for the detailed dynamic nonlinear processes involved in disease transmission, also sometimes referred to as “process-based”), have been employed to assess the likely impact of anthropogenic climate change on malaria transmission dynamics and control. These models have reached divergent conclusions, with some predicting a large expansion in the continental land area suitable for transmission (Martens et al. 1999; Caminade et al. 2014; Tanser et al. 2003) and in the number of people at risk of malaria (Martens et al. 1999; Patz et al. 1996; Pascual et al. 2006), while others predict only modest poleward (and altitudinal) shifts in the burden of disease, with little net effect (Gething et al. 2010; Rogers and Randolph 2000; Hay et al. 2002), and the issue remains unresolved thus far. The goal in this paper is not to ultimately resolve this issue (laudable as it is), but to attempt to lay a foundation to aid such a resolution.

Malaria was one of the first human diseases to be subject to mathematical inquiry. Sir Ronald Ross, who first elucidated how Plasmodia spread via the intermediary mosquito (Cox 2010), proposed a mechanistic transmission model including the human host and the mosquito in the early 1900s, although it did not address the all-important mosquito lifecycle (beyond infection of a constant population) (Smith et al. 2012). Following in Ross’s footsteps, the highly influential malariologist George Macdonald reformulated the basic model in the early 1950s (Macdonald 1952, 1956a, b, 1957) (presented in detail in Sect. 4.2), and derived an expression for the basic reproduction number, ${\mathcal R}_0$, defined as the average number of secondary cases a single initial case will generate in a completely susceptible (uninfected and non-immune) population.¹ Macdonald showed that ${\mathcal R}_0$ is most sensitive to changes in adult mosquito survival probability, thereby providing a theoretical rationale for insecticide spraying as the foundation of malaria eradication efforts during the 1960s (Macdonald 1956b; Nájera et al. 2011). The “Ross–Macdonald” model has been extremely influential: Reiner et al. (2013), in a systematic review of 388 models of mosquito-borne pathogens found in the literature between 1970 and 2010, determined most to be similar to the Ross–Macdonald framework, and Macdonald’s expression for ${\mathcal R}_0$ has also been used in many climate-focused (or climate-driven) mechanistic modeling studies.

Moving back to the effect of climate change, there has been significant controversy as to its likely effect on human malaria (and other diseases) and this has been informed by modeling studies. Several works in the 1990s by Martens, Lindsay and colleagues (Martens et al. 1995a, b; Lindsay and Martens 1998; Martens et al. 1997, 1999), using Macdonald’s ${\mathcal R}_0$, and drawing from several sources quantifying how parasite and vector lifecycle parameters might change with temperature, predicted that a significantly expanded area of the globe could become vulnerable to epidemic malaria under climate change. However, these conclusions did not go unchallenged.

Rogers and Randolph (2000) were critical of these process-based methods, and instead employed a statistical method, whereby they inferred climatologic limits to malaria based on temperature, rainfall, and saturation vapor pressure, and current malaria distribution, and then projected how malaria suitability would change in the future, under global climate model (or general circulation model, GCM) projections, finding a decrease in some mainly equatorial areas and a modest poleward increase, with little net change overall; this agrees in principle with a historical time-series analysis by Small et al. (2003). Other authors have similarly argued that climate change is more likely to induce a geographic shift in the burden of disease, with little net increase (Lafferty 2009). Even if true, Pascual and Bouma (2009) have pointed out that a geographically balanced shift does not equal a population-balanced shift: highland regions in eastern Africa, most likely to become vulnerable to malaria under a warming climate, are also far more populous than nearby lowland areas that could see a decline in malaria burden.

Gething et al. (2010) also argued, essentially, that because the global burden of malarial disease decreased dramatically from 1900 to 2007 while global mean air temperatures increased (the average temperature increase from the 1850–1900 period to the 1986–2005 period was 0.61 $^{\circ }\hbox {C}$ (IPCC 2014)), then non-climactic factors must be of vastly overriding importance and that climate change will affect malaria but little in the future. They bolstered this argument by estimating how ${\mathcal R}_0$ must have changed overall since 1900 and in response to different interventions, based on a Ross–Macdonald-style model for ${\mathcal R}_0$, and concluded that projected mean increases in ${\mathcal R}_0$ under future warming (Martens et al. 1997, 1999; Lindsay and Martens 1998) are one to two orders of magnitude smaller (and, thus, likely trivial). Agricultural practices and land-use, in particular, have been considered a human factor of greater importance than climate change (Lafferty 2009).

A variety of newer process-based models for the transmission cycle were developed in the decade after Martens, varying in their basic construction and hypotheses for the effects of rainfall and temperature on both vector and parasite (see, for instance, (Hoshen and Morse 2004; Bomblies et al. 2008; Parham and Michael 2010; Alonso et al. 2011; Ermert et al. 2011a, b; Parham et al. 2012; White et al. 2011)), but generally concluded that increasing temperatures favor malaria transmission. For example, Parham and Michael (2010) concluded in 2010 that transmission is optimized in the 32–33 $^{\circ }\hbox {C}$ temperature range. Caminade et al. (2014) published projections for the population at risk of malaria using five malaria models from this period, suggesting a net increase in the global population at risk of malaria, but with high uncertainty.

Mordecai and colleagues (Mordecai et al. 2013), in an influential paper published in 2013, used a set of unimodal functions (i.e., hump-shaped) for the relationship between temperature and vector parameters (such as larval development rate, larval survival, adult survival, biting rate, fecundity, and vector competence), as well as the parasite development rate, in contrast to many prior works, which had used monotonic relationships for some or all of these (temperature-dependent) parameters. Using a newer expression for ${\mathcal R}_0$, based on the model by Parham and Michael (2010), Mordecai and colleagues concluded that malaria transmission is optimized at a significantly lower temperature range, 25–28 $^{\circ }\hbox {C}$, and found this to better match field measurements of the entomological inoculation rate (EIR) (the infectious biting rate).

Subsequently, Ryan et al. (2015b) used the Mordecai et al. (2013) thermal-response curves to develop a series of maps for malaria transmission potential across Africa from the year 2000 to 2080 under a mid-range emissions scenario (SRES A1B). Broadly speaking, this work predicted a modest increase in the total land area at risk for any malaria transmission, while the net area suitable for intense, year-round transmission would decrease (especially in western coastal Africa). Furthermore, these authors predicted increased malaria potential in the cooler southern and eastern regions of Africa, but a decrease in the hotter western and central African regions (especially the Democratic Republic of Congo) by 2080, and an southeasterly shift over time in the populations most at risk of malaria, with notable increases in the Lake Victoria region (near the Kenyan highlands) and eastern highland Madagascar. This work is especially laudable in its nuanced approach to malaria transmission potential, differentiating between year-round and seasonal potential, and consideration of the populations, not just geographic areas, at risk.

Despite its virtues, the work of Ryan et al. (2015b) did not explicitly consider rainfall or hydrodynamics, but applied a mask that limited transmission only to those regions with enough vegetation to be considered wet enough to support anopheline habitat. Earlier (process-based) malaria potential maps based on temperature, e.g. that of Craig et al. (1999), somewhat similarly restricted transmission to areas with grossly sufficient rainfall. Indeed, most of the works reviewed thus far have focused primarily on ambient temperature as an explanatory variable, with rainfall often a secondary, and variously modeled, factor. Given the absolute necessity of appropriate aquatic habitat to the vector lifecycle, hydrodynamics and habitat modeling at both the regional and micro-scale represent a relatively (but not entirely) neglected factor. A variety of relatively simple relations between rainfall and immature mosquito survival and carrying capacity have been employed (Yé et al. 2009; White et al. 2011; Hoshen and Morse 2004), while several more complex efforts (Paaijmans et al. 2008a, b; Parham et al. 2012; Asare et al. 2016a, b) have physically modeled the heat and water balance within Anopheles microhabitats, as reviewed in Sects. 5.2.5 and 5.2.6. Several authors have additionally modeled regional hydrodynamics, e.g. (Bomblies et al. 2009; Bomblies 2012; Tompkins and Ermert 2013; Asare et al. 2016c). Of especial note, Bomblies and colleagues have considered detailed hydrodynamics at the village scale (Bomblies et al. 2008, 2009; Bomblies 2012), and concluded that such detailed modeling is necessary to explain both interseasonal variation (Bomblies 2012) and intervillage variation in vector abundance (Bomblies et al. 2009), and this modeling formed the basis for a recent comprehensive study suggesting little effect of climate change on malaria incidence in western Africa (Yamana et al. 2016).

While much controversy has centered on the appropriate functions relating vector and parasite parameters to temperature (and secondarily, to rainfall) and how variations in these drive climate-related predictions, more basic modeling choices also affect model predictions. In particular, the population biology of the Anopheles vectors is crucial to understanding many aspects of the disease, as well as assessing control strategies and projecting future outcomes. Malaria models that do not incorporate the dynamics of the juvenile stages of the mosquito are known to give results that do not generally match observed epidemiology (Okuneye and Gumel 2017; Beck-Johnson et al. 2013), and the vector lifecycle per se is the focus of several models (Beck-Johnson et al. 2013), most recently by Abdelrazec and Gumel (2017), who studied the effect of both temperature and rainfall on the population biology of mosquitoes. Another fundamental issue is that most vector and parasite lifecycle process times (e.g., larval development time) are non-exponentially distributed, yet most differential equations-based disease transmission models implicitly assume exponentially-distributed waiting times, an assumption found to affect model dynamics unfavorably by Christiansen-Jucht et al. (2015) and Lunde et al. (2013b).

Addressing this deeper problem of model construction, Gumel and colleagues (Agusto et al. 2015; Okuneye and Gumel 2017), have recently developed and analyzed several complex weather-driven mechanistic models that extend the prior studies by incorporating a broader array of biological, ecological and epidemiological factors, such as the dynamics of immature mosquitoes, host age-structure (Okuneye and Gumel 2017) and host immunity-boosting due to repeated exposure to malaria infection (Agusto et al. 2015). In particular, Agusto et al. (2015), adopting the thermal-response functions of Mordecai et al. (2013), and using a 14-dimensional mechanistic model and weather data for numerous locations within Africa, predicted that malaria infection generally increases in the 16–28 $^{\circ }\hbox {C}$ range, but decreases beginning at temperature values between 25 and 28 $^{\circ }\hbox {C}$, depending on the African region (these results are comparable to those of Mordecai et al. (2013), but more nuanced). Yamana et al. (2013) also extended a prior agent-based model by Bomblies et al. (2008) to include partial immunity induced by repeated infection, and predicted that immunity can damp both the spatial and temporal variation in clinical disease in response to environmental variability (Yamana et al. 2013, 2017). It should be emphasized that many prior weather-driven malaria modeling studies do not include immunity (or use only very simple representations of immunity), even though it is known that the unique malaria immune response is fundamental to malaria epidemiology and pathogenesis, and is itself the focus of a long modeling tradition; see, for instance, (Dietz et al. 1974; Aron 1988; Gupta and Day 1994; Gupta et al. 1999a, b; Filipe et al. 2007; Griffin et al. 2010, 2015) and Sect. 7.1.

Another recent effort is that of Okuneye and Gumel (2017), who additionally incorporated age-structure (as stated earlier, age-structure is crucially important because children under the age of five suffer the majority of the malaria burden in endemic areas) into a mechanistic temperature- and rainfall-dependent model, finding transmission to be maximized in the 21–25 $^{\circ }\hbox {C}$ temperature and 95–125 mm rainfall ranges in the Kwa-Zulu Natal province of South Africa.

Yet another basic issue that must be mentioned is that of diurnal temperature variation, and the (time-varying) disparity between ambient air and water temperature. Paaijmans et al. (2010) demonstrated empirically that the magnitude of temperature fluctuation affects Anopheles development and survival in a manner not captured by mean temperature alone. Average diurnal temperature range varies on a continental scale (Paaijmans et al. 2010), and this therefore may be an under-appreciated parameter in malaria potential projections. Diurnal temperature variations have not been considered in most models, although there are some recent exceptions, e.g. (Agusto et al. 2015; Beck-Johnson et al. 2017). Furthermore, the water temperature in immature mosquito habitats generally differs from ambient air; this disparity may be captured by physical hydrodynamic modeling, although a simple linear offset is sometimes assumed (Agusto et al. 2015). Finally, adult anophelines are also exposed to multiple microenvironments with varying temperatures, and often prefer to feed and/or rest indoors, where temperatures are typically warmer on average, but also less variable than out-of-doors (Afrane et al. 2005; Blanford et al. 2013; Singh et al. 2016).

While many malaria modeling studies have focused on the global scale (i.e., the potential global malaria range due to climate change), studies more limited in scale may provide better insight (Pascual and Bouma 2009; Alonso et al. 2011). In particular, a model region is the highlands of East Africa, where malaria burden was previously rare but has become more common since the 1970s; this increase may be at least partially attributable to global warming (Pascual and Bouma 2009). Human activity in the Kenyan highlands is recapitulating, in some sense, the early social and climatic changes that first gave birth to P. falciparum some 10,000 years ago. Temperatures are increasing (Pascual et al. 2006), the rain forests have recently been mostly cleared for crops, cattle grazing, logging, and housing construction (Minakawa et al. 1999), and the region is subject to intense population growth and human migration. Several researchers have made this area their focus (e.g., Githeko and Ndegwa 2001; Hay et al. 2002; Zhou et al. 2004; Pascual et al. 2006, 2008; Chaves and Koenraadt 2010; Alonso et al. 2011; Snow et al. 2015), and we suggest that a more limited geographic scope of study may better elucidate the competing effects of treatment, land use, migration, and climate on malaria. Also of note, malaria is highly endemic in hotter western Africa, an area which is also the focus of several studies, and the effect of climate change in this region could, conversely, be to slightly reduce malaria potential (Ryan et al. 2015b; Yamana et al. 2016).

In summary, although there is general agreement that climate change will increase the potential for malaria transmission at more northerly and southerly latitudes (and at higher altitudes), it is unclear if this represents a shift in malaria distribution with little net increase (or even decrease) in malaria burden, or an expansion in burden. The most likely scenario may be a hybrid result, with net expansion in malaria range, but shifts in the intensity of transmission within that range, especially towards southern and eastern Africa and highland areas. Further, the magnitude of the climate effect, and how it compares to other anthropogenic and abiotic factors, remains uncertain.

Malaria is a complex disease, with a complex history, and the controversy just outlined cannot be fully addressed without a broad background. The goal of this paper is to provide the reader with at least some of the requisite background needed for effective modeling of the disease dynamics, and to provide sufficient resources to help the reader in beginning their own investigations. Finally, it should be noted that mathematical models can, broadly speaking, be divided into the classifications of non-parametric and parametric (with parametric models also referred to as “process-based” or “mechanistic”), where the former attempt to inferentially draw conclusions directly from (usually time-series) data without positing any particular mechanistic system, while the latter posit some particular hypothesis for a system’s workings (expressed mathematically). In this paper, our focus shall be on the latter.

This paper is organized as follows. We begin with an overview of the malaria lifecycles, immunology, and epidemiologic principles to establish a basis for later sections. To properly appreciate the role of mathematical modeling in providing deeper qualitative and quantitative insight on the transmission dynamics and control of malaria, some familiarity with the historical development of modeling frameworks and concepts is invaluable. To this end, we first present a historical overview of the disease in general, and move on to quantitative malariology through the early twentieth century, focusing on the early but deeply influential work of Ross and Macdonald. We also touch on some important later extensions by authors including Garrett-Jones, Dietz, and Molineaux. We then shift focus to climate, beginning with an extensive discussion of the anopheline and parasite lifecycles and their relation to weather (mainly temperature and rainfall), since these are fundamental to any predictions we care to make. In Sect. 6, we subsequently present a partial genealogy of recent mathematical works addressing weather and malaria transmission, and close with a brief discussion of multi-patch meta-population modeling, which may be of especial importance in understanding the spread of malaria between lowland and highland regions of Kenya. Finally, we briefly discuss other aspects of the disease that are pertinent to a fully comprehensive quantitative modeling framework (such as malaria immunity).

2 Introduction to malaria lifecycles, immunology and clinical disease

2.1 Parasite lifecycle

Plasmodium spp., the causative agent in malaria, are sexually-reproducing eukaryotic protozoans that undergo a complex lifecycle that requires switching between evolutionarily-distant vertebrate and invertebrate dipterian hosts. The basic evolutionary logic follows. Pre-Plasmodium parasites likely evolved from free-living sexual protozoans to live extracellularly in the midgut of aquatic invertebrates (Carter and Mendis 2002). Proliferative potential was then increased with the evolution of a second parasitic intracellular asexual reproductive stage, known as schizogony, by which a single cell may produce vast numbers of daughter cells, or spores. A minority of these daughter spores differentiate into male and female forms, which then recombine in a form of extracellular sexual reproduction known as sporogony (Antinori et al. 2012). The Plasmodia’s evolutionary innovation is to spatially separate the schizogonic cycle and sporogonic cycle into two separate hosts, with sporogony occurring in the mosquito.

Let us consider the particulars of human Plasmodia, where schizogony (asexual clonal expansion of many daughter spores) occurs in the human host, and in two phases: first in liver hepatocytes and then within red blood cells (RBCs, or erythrocytes). Sporogony then occurs in the mosquito midgut following a blood meal, to ultimately yield parasitic forms infectious to humans (Antinori et al. 2012). We may consider the cycle to begin with the bite of an infectious mosquito, who probes the dermis and injects saliva containing no more than 10–100 highly motile asexual sporozoites (Antinori et al. 2012). Sporozoites penetrate into blood vessels within minutes, travel to the liver and establish infection in hepatocytes within 30 min of biting (Guilbride et al. 2012). While the skin has traditionally been thought of as a passive waypoint in the infection cycle, more recent data indicates that a small number of sporozoites remaining in skin may exploit the inherently immunoregulatory nature of this environment to suppress anti-Plasmodium immunity and induce tolerance (Crompton et al. 2014), with important implications for vaccine development (Guilbride et al. 2012).

Shifting focus, sporozoites within hepatocytes initiate the first round of shizogony, so-called “pre-erythrocyte” shizogony, proliferating asexually to produce, in the case of P. falciparum, up to 30,000–40,000 asexual merozoites (Antinori et al. 2012; Crompton et al. 2014) contained within a “tissue schizont”. Once mature, the tissue schizont, along with the parent hepatocyte, ruptures to spill the merozoites into the bloodstream, where they actively infect red blood cells, initiating the erythrocyte cycle of schizogony (Antinori et al. 2012), whereby merozoites expand, via several intermediate stages, within the erythrocyte and rupture it every 24–72 h (48 h for P. falciparum), freeing more merozoites to repeat the cycle (Antinori et al. 2012). It should be noted that while this is the end of the hepatic stage for P. falciparum, P. vivax and P. ovale have a dormant liver form known as the hypnozoite (Greek “sleeping animal”) that can cause reinfection years later (Carter and Mendis 2002).

Erythrocytes, lacking a nucleus and most typical eukaryotic organelles, are essentially masses of hemoglobin, an iron-containing oxygen-carrying molecule, wrapped in plasma membrane and suited only for passive O$_2$ and CO$_2$ transport. Plasmodia, on the other hand, are “fully realized” eukaryotes, that hijack completely the erythrocytes they invade (Tilley et al. 2011). An invading merozoite passes through an immature “ring” stage to become a trophozoite, a feeding form that consumes 70% of the erythrocyte hemoglobin, converting it to the toxic byproduct hematin, which is then detoxified to hemozoin (Baton and Ranford-Cartwright 2005; Tilley et al. 2011). Notably, quinine antimalarial drugs act by preventing the detoxification of hematin (Tilley et al. 2011), and artemisinin, the most effective antimalarial, is also likely involved in hemoglobin digestion (Tilley et al. 2011). The nourished trophozoite then becomes a “blood schizont,” dividing asexually into 6–36 (20 on average) daughter merozoites that are released with the host cell’s rupture (Tilley et al. 2011).

The erythrocyte cycle can continue essentially indefinitely, but it is ultimately a reproductive dead end: the merozoite must die with the man. To escape the human host and live on, a subset of blood schizonts commit their merozoite offspring to becoming gametocytes, sexually differentiated male and female parasite forms; all progeny of a schizont become either male, female, or asexual (the most typical fate). Sexually committed merozoites proceed as others, by invading an erythrocyte to become a feeding trophozoite, but then form either a single macrogametocyte (female) or single microgametocyte (male) (Baton and Ranford-Cartwright 2005). Committing to gametocytogenesis is risky, for terminally differentiated gametocytes cannot reproduce further in the blood of man. Gametocytogenesis in Plasmodia does not occur after a fixed number of erythrocyte cycles, as it does in some related parasites, and the decision to commit to gametocytogenesis remains poorly understood (Baton and Ranford-Cartwright 2005).

When taken up in a blood meal, the gametocytes rapidly initiate the sporogonic cycle. Upon arrival at the mosquito midgut the macrogametocyte dissolves its erythrocyte host within minutes and becomes spherical and immotile. The microgametocyte, on the other hand, undergoes the dramatic process of exflagellation, whereby eight daughter genomes are produced that attach to long writhing flagella, and break free to become highly motile wrigglers that find and fuse with a macrogametocyte to form the zygote (Baton and Ranford-Cartwright 2005; Antinori et al. 2012). The zygote in turn transforms into a banana-shaped ookinete (Greek “moving egg”), which penetrates both through the peritrophic matrix, a chitinous matrix extruded by the mosquito gut to sequester the blood meal, and then through the epithelial cells lining the gut wall. Next, it transforms into an oocyst, producing thousands of daughter sporozoites by nuclear division (Antinori et al. 2012). Eventually, the mature oocyst ruptures into the mosquito’s hemocoelic cavity (Baton and Ranford-Cartwright 2005), and sporozoites travel through the hemolymph to infect the salivary glands where, after about one day, they are reprogrammed to be highly infectious to humans (Antinori et al. 2012), and the cycle can begin again. The cycle is depicted in its entirety in Fig. 2.

The complex within-host dynamics of human Plasmodium infection, how these affect the efficacy of treatment and control measures, and their interaction with the immune response, have been the focus of multiple modeling works, e.g. (Teboh-Ewungkem et al. 2010; Li et al. 2011; Gurarie et al. 2012; Eckhoff 2012; Demasse and Ducrot 2013; Childs and Buckee 2015; Childs and Prosper 2017; Tabo et al. 2017). However, to our knowledge no climate-focused models have focused deeply upon these within-host dynamics, although it is likely that such work is needed to fully elucidate how climate change might affect malaria epidemiology and control efforts in the future (see also Sect. 7.3).

Fig. 2

The Plasmodium lifecycle. The right side depicts schizogony in man, where sporozoites from an infectious bite invade hepatocytes in the liver, undergo a round of replication, and then enter, as merozoites, into the erythrocyte cycle in blood. A minority of blood trophozoites differentiate to male and female gametocytes that are taken up by a biting mosquito to initiate sporogony, as depicted on the left side, whereby ookinetes penetrate the peritrophic matrix and gut epithelium to form oocysts, eventually rupturing to yield sporozoites that travel to the salivary glands. Further details are given in the text

2.2 Vector characteristics and lifecycle

Malaria is transmitted by adult female Anopheles mosquitoes, yet of the more than 450 known anopheline species, only about 60 can serve as actual vectors (Cohuet et al. 2010), with 41 considered major vectors (Sinka et al. 2010), and most of these are rather inefficient at transmitting the disease. To effectively transmit disease, the mosquito must be susceptible to infection (many are completely refractory to Plasmodium), must habitually bite man (many mosquitoes strongly prefer other animals), and must live long enough for the sporogenic cycle to reach completion (Cohuet et al. 2010).

In Africa, three anopheline species are preeminent, namely A. gambiae, A. arabiensis, and A. funestus, with A. gambiae likely the single most important species (Sinka et al. 2010), and the focus of most modeling studies. A point of terminology to avoid confusion in the literature is in order here: the A. gambiae complex is a collection of seven morphologically indistinguishable species later recognized to be distinct, and includes A. gambiae sensu stricto (Latin “in the strict sense”) which is the species referred to by the unqualified term A. gambiae, and A. arabiensis (Sinka et al. 2010). A. gambiae sensu lato (Latin “in the general sense”) refers to the species complex. The existence of multiple distinct species, including the important vector A. arabiensis, within the A. gambiae complex clearly complicates matters, from both a modeling and malaria control standpoint, and these vectors vary, for example, in their susceptibility to insecticide-treated bednets (Kitau et al. 2012).

Briefly, the lifecycle of the Anopheles mosquito, while simpler than that of its Plasmodium parasite, is not trivial, with mosquitoes passing through three immature, aquatic stages (egg, larva, pupa), and a final adult stage. The adult female mosquito lifecycle is centered around the gonotrophic cycle: the taking of a blood meal to fuel egg development, which takes several days and is highly temperature dependent (with higher temperatures decreasing the time for larval development), followed by oviposition of eggs in a suitable aquatic habitat, only to repeat until inevitable death (Detinova 1962). Typical blood meal size is 2–3 $\upmu$L, and A. gambaie may lay anywhere between about 10 and 150 eggs per gonotrophic cycle (this is the “fecundity”) (Takken et al. 1998; Afrane et al. 2005), but more typically between about 40 and 85 under field conditions (Afrane et al. 2005). Most eggs hatch within 2–3 days (Yaro et al. 2006), but time to hatching is modestly temperature dependent (Bayoh and Lindsay 2003). Eggs hatch to become larvae and actively feed upon algae and bacteria, growing through four moltings, and are thus divided into four stages known as instars (conceptually, first- and second-instars are lumped as “early,” with third- and fourth- “late-instars”); Anopheles larvae also lie parallel to the water surface to obtain oxygen. Finally, fourth instar larvae become nonfeeding pupae that undergo metamorphosis to adult mosquitoes. Immature stage development rate and survival are both strongly temperature-dependent (Bayoh and Lindsay 2003), and the complete lifecycle is given in schema in Fig. 3.

Fig. 3

Anopheles mosquito lifecycle. Immature mosquitoes pass through aquatic egg, larvae, and pupae stages, with the actively feeding larvae divided into four instar stages. Adult female mosquitoes pass through the gonotrophic cycle, by which bloodmeals nourish the development of new eggs, with further details in the text

Anophelines have varying habitat preferences, and are widely adapted to different environmental niches (Sinka et al. 2010), but the A. gambiae complex tends, unsurprisingly, to prefer conditions associated with anthropogenic alteration of the environment. Specifically, A. gambiae and A. arabiensis larvae prefer small, temporary, sunlit pools, with little vegetation (Minakawa et al. 1999, 2004), the kind created by deforestation, construction, and livestock, e.g. hoofprints. These pools are warmer, support more algae (the major larval food source), and have fewer predators than natural water bodies (Minakawa et al. 1999). A series of studies by Afrane and colleagues (Afrane et al. 2005, 2007, 2008) confirm that deforestation in Kenyan highlands creates habitat that strongly supports A. gambiae proliferation. A. funestus, the other major African vector, is also aided by deforestation, but tends to prefer larger permanent or semipermanent habitats with established vegetation (Minakawa et al. 2005).

Anophelines can further be characterized along several axes that affect transmission potential and the efficacy of various control efforts (Sinka et al. 2010): (1) anthropophilia versus zoophilia, or the preference for taking blood meals from humans or non-human animals, respectively, with the anthropophilic index defined as the fraction of blood meals taken from man, (2) endophagic versus exophagic, referring to a preference for feeding indoors or out-of-doors, respectively, and (3) endophilic versus exophilic, meaning the favored location for resting between blood meals (this may differ pre-feeding and post-feeding). The highly efficient malaria vectors, such as A. gambiae, tend to be highly anthropophilic with anthropophilic indices approaching one, but as reviewed by Sinka et al. (2010), even these vectors are likely very opportunistic, and apparent anthropophilia may simply be a (partial) result of host availability; preferences along the other axes may also have been overstated in past studies, and anophelines are quite adaptable in general (Sinka et al. 2010).

Briefly, A. gambiae feeds late at night, has typically been reported as endophagic and endophilic, but this likely varies, and Odiere et al. (2007), working in western Kenya, found no preference. The closely related A. arabiensis also feeds at night, and may show an exophagic and exophilic preference in comparison to A. gambiae (Sinka et al. 2010). Behaviorally, the adult A. funestus is extremely similar to A. gambiae (Sinka et al. 2010).

Finally, not only do the innate characteristics of certain anophelines favor malaria spread, but there is even some evidence that mosquito behavior may also be modulated by Plasmodium infection to enhance transmission, a notion termed the “manipulation hypothesis” by Cator et al. (2012). When carrying the infectious sporozoite parasite stage, various Anopheles may take more frequent bloodmeals with more probing attempts per meal, may be more likely to feed from multiple hosts, and bloodmeal volume may be smaller. In contrast, when burdened by the non-infectious oocyst stage, mosquitoes seem less attracted to hosts and less persistent in bloodmeal attempts, a behavioral response that could decrease pre-infectious mortality by avoiding risky biting attempts (Cator et al. 2012; Nguyen et al. 2017), and the overall effect of such manipulations on malaria transmission could be quite significant, as suggested by mathematical analysis by Cator et al. (2014). However, most evidence for such manipulation comes from lab studies using a variety of vector-host combinations (Cator et al. 2012), and it is also unclear if such behavioral changes represent specific parasitic manipulations or more generic responses to infection (Cator et al. 2013). Moreover, several recent studies using field isolates of P. falciparum and anthropophilic Anopheles found no evidence that infection altered host-seeking behavior (Vantaux et al. 2015; Nguyen et al. 2017).

2.3 Immunity and clinical disease

In areas of intense P. falciparum transmission, young children are exposed to hundreds of infectious bites per year, and yet, unlike many viral diseases where a single exposure can be sufficient to imbue robust, lifelong immunity, immunity to malaria is gained only slowly and incompletely over the course of years (Crompton et al. 2014). Characteristically, children under the age of five are susceptible to the most severe, life-threatening forms of the disease, such as severe malarial anemia and cerebral malaria, the disease transitions to an uncomplicated febrile disease through adolescence, and by adulthood it only rarely manifests clinically, with asymptomatic disease common (Crompton et al. 2014) (and a possible Plasmodium reservoir complicating eradication efforts). This hard-won immunity is short-lived: when adults from endemic areas move, they become vulnerable to severe disease within a few years, although they may retain protection against the worst manifestations of disease (Filipe et al. 2007). This dynamic is especially salient to malaria elimination efforts, and there is very real danger when the disease is eliminated locally but may be reintroduced to now non-immune populations (Webb 2014; Snow 2015), and even control measures, such as bednets or intermittent preventive (drug) therapy, while initially beneficial, have the potential to increase disease burden later in time, as they induce a decrease in population-level immunity (Ghani et al. 2009).

Thus, it is clear that there is a distinct disparity between clinical immunity against P. falciparum malaria (protection against clinical disease and severe symptoms) and infectious immunity (protection against infection, per se, by blood-stage parasites). Immunity to the most severe forms of disease may be also differ fundamentally from immunity to uncomplicated disease, with perhaps only several infections (and possibly as few as one) needed to confer long-lasting protection (Gupta et al. 1999a, b). The pathogenesis of clinical disease is primarily related to (1) sequestration of parasitized erythrocytes in organs such as the brain (this sequestration allows parasites to avoid the spleen, where they could be destroyed by macrophages), and (2) the systemic inflammatory response (Crompton et al. 2014). With respect to the former, the P. falciparum erythrocyte membrane protein-1 (PfEMP1), expressed on the surface of infected cells, facilitates sequestration. It is also encoded on the var gene, of which there are about 60 distinct versions, each clonally expressed and encoding an antigenically distinct PfEMP1. This antigenic variation, and the extreme genetic diversity of P. falciparum in general, help to explain why effective immunity requires so many exposures (Crompton et al. 2014).

It is worth noting that all actively clinical disease takes place during the erythrocyte stage of infection, with the skin and hepatocyte stages clinically silent. This may be at least partly related to the very different orders of magnitude involved at the different stages. Generally, fewer than 100 sporozoites infect the skin, and only tens of hepatocytes are infected. These numbers may simply be too small to initiate immunity, or, they may even induce immune tolerance, especially in the skin (Guilbride et al. 2012). In severe infections, on the other hand, total body trophozoite burden may number in the hundreds of billions (Trape et al. 1994).

2.4 Epidemiologic classification

P. falciparum transmission intensity in endemic zones varies across orders of magnitude, from one infectious bite per person per year, to more than one per day in many holoendemic areas (Rodriguez-Barraquer et al. 2016), and partly as a consequence of its unique immunology, different malaria transmission intensities give differing age-distributions of parasitemia, clinical disease, and serious disease (Aron 1988; Snow 2015). It must be emphasized that endemic and epidemic malaria are very different beasts (Snow 2015; Hay et al. 2008): endemic (from Greek meaning “in the people”) disease is constantly present in a population, whereas an epidemic (Greek “upon the people”) is a temporary disease flare out of proportion to the past. Populations living with endemic malaria have varying degrees of immunity, but epidemic malaria can be calamitous when it tears through previously unexposed groups, or more perniciously, groups transiently protected by malaria control programs that lapse, leaving the people newly vulnerable after the waning of prior immunity (Snow 2015; Webb 2014).

The most common classification for endemicity is now based upon the fraction of the population that has parasites detectable in their peripheral blood, the so-called “parasite rate” (PR), and furthermore uses the parasite rate in the 2–10 year age group, PfPR$_{2-10}$, with zones classified as holoendemic (PfPR$_{2-10}$ > 75%), hyperendemic (PfPR$_{2-10}$ 50–75%), mesoendemic (PfPR$_{10}$ 10–50%), and hypoendemic (PfPR$_{10}$ 1–10%) (Snow 2015). In the hypoendemic and extreme hypoendemic (< 1%) range, transmission becomes unstable (Snow 2015), and populations with very low burdens of malaria are vulnerable to epidemics of severe disease. Indeed, the venerable Macdonald considered the stable/unstable classification axis to be the more legitimate on a fundamental basis (Snow 2015).

The age-distribution of clinical disease varies across endemic zones. In holoendemic zones, most severe disease occurs in the first few years of life, rapidly tapering off by adolescence (Aron 1988; Gupta and Day 1994; Filipe et al. 2007; Crompton et al. 2014; Snow 2015), with the burden of severe disease dropping in absolute terms and shifting towards older age groups as the level of endemicity decreases, as demonstrated in Fig. 4. In holoendemic areas, the PR peaks later than does clinical disease (Trape et al. 1994; Rodriguez-Barraquer et al. 2016), and remains relatively high even into middle and old age, when clinical disease is rare. However, although the PR remains high, the parasite burden continues to decline with age (Trape et al. 1994), as also shown in Fig. 4. These observations have motivated many mathematical models attempting to elucidate the dynamics of immunity acquisition.

Fig. 4

The left panel gives the qualitative shape of severe disease incidence by age through adolescence under holo-, hyper-, meso-, and hypoendemic transmission conditions, based on Aron (1988), Snow (2015). The right panel shows overall parasite rate in the holoendemic village of Dielmo, Senegal (Trape et al. 1994), subdivided by the actual density of trophozoites in the peripheral blood. While PR in the youngest is only about twice that of those over 40, children under four suffered clinical malaria attacks at a rate 40-fold higher (Trape et al. 1994). Note that the broad plateau in PR from roughly age two to 15 in the face of dramatically falling serious disease incidence has been observed elsewhere (Gupta and Day 1994)

3 General historical background

3.1 Overview

The history of malaria, and its emergence as a major human pathogen over the last several 10,000 years, is intimately linked to the evolution of human agricultural civilization and the profound changes in both human populations and the environment that this engendered. This was directly coupled to global climate, as climate change following the end of the last ice age and the onset of the holocene era was fundamental to agriculture, and also allowed the wider spread of mosquito vectors in a warmer world (Carter and Mendis 2002). For the interested reader, scholarly histories of the disease include those by Webb (2014) and Packard (2007).

The clearing of forests for agriculture creates myriad microenvironments for anopheline mosquitoes, and concentrated human settlements are capable of supporting the virulent P. falciparum, which only emerged within the last 10,000 years, while a warmer climate supports its mosquito vectors (Webb 2014; Packard 2007). Malaria has also helped shape human biological evolution: in pre-agricultural Africa Duffy antigen (an erythrocyte membrane chemokine) negativity spread through the African heart to confer complete resistance to P. vivax, at no apparent cost, while the more recent advent of P. falciparum selected for a variety of far less benign genetic anemias, the best known being the sickle cell trait, which protects against severe disease in the heterozygous form, but causes crippling sickle cell disease in homozygotes (Carter and Mendis 2002). Following its earlier evolution in Africa, malaria, especially P. vivax, escaped that continent and into much of the rest of the world, its spread strongly associated with agricultural expansion and population movements (Packard 2007).

It was not until the end of the nineteenth century that the microbiological basis of the disease was discerned. This coincided with the onset of the colonial era, or “Scramble for Africa” spanning roughly 1879 through 1914, and during which various European powers conquered and carved up the African continent (Webb 2014). Thus, early “scientific” malaria control efforts in Africa were inescapably linked to colonial medicine, a primary focus of which was protecting Europeans and preserving the productivity of subservient African laborers, with less regard for the general African populace (Webb 2014). Lasting from 1955–1969, the World Health Organization’s Malaria Eradication Programme saw significant mixed successes, but ultimately failed to eliminate the disease. In Africa, widespread chloroquine treatment during the 1970s was a primary cause of lowering malaria burden, but the spread of chloroquine resistance across the continent in subsequent years, the HIV/AIDS epidemic, agricultural expansion, and devastating wars among many newly independent African states fueled a malaria resurgence. New global efforts since about 2000, largely centered on insecticide-treated bednets (ITNs) and treatment with the new artimisinin compounds have seen significant success (Bhatt et al. 2015), but it remains to be seen whether these gains will continue or even be maintained (Webb 2014).

In the following sections, we discuss more extensively the early origins of malaria, its link to agriculture and human activity, and then review in greater depth the era since the late nineteenth century. Our focus in the latter is on tropical Africa and P. falciparum, and a working theme is that P. falciparum differs qualitatively from the other human Plasmodia, representing a unique burden on African populations.

3.2 Origins and evolution

The Plasmodia are ancient parasites belonging to the order haemosporidia—single-celled parasites which alternate between a wide variety of vertebrate hosts and blood-sucking arthropods—and mammal-specific Plasmodia have coexisted with mammals for much if not all their evolutionary history, with one recent estimate dating their origin between 64 and 120 million years ago (Silva et al. 2015). Haemosporidia burdened animals even earlier, likely since almost the first appearance of Diptera insects (flies and mosquitoes) 150–200 million years ago (Carter and Mendis 2002).

Early studies found P. falciparum and a very closely related chimpanzee Plasmodium, P. reichenowi, to differ substantially in morphology and lifecycle from P. malariae, P. ovale, and P. vivax, and hence the former were categorized as a separate subgenus, Laverania (Loy et al. 2017). Later molecular studies confirmed that the Laverania diverged from the other mammalian Plasmodia on the order of 100 million years ago (Carter and Mendis 2002; Silva et al. 2015).

Moving forward in time, the evolutionary origins of modern human Plasmodia within the last 100,000 years, mainly P. vivax and P. falciparum, have been of some controversy, but the weight of the evidence supports, in our view, an out-of-Africa origin for all modern human malaria (see Loy et al. (2017) for a recent review). Under pre-agricultural conditions, scattered mobile populations of low density were unlikely to support intense transmission rates, and the overall malaria burden was probably low. Under these poor transmission conditions, P. malariae, which can cause a chronic low-grade infection lasting decades, and P. vivax and P. ovale, both of which have a dormant liver stage that can lead to reinfection and transmission years after initial infection, are much more competitive than the highly virulent and short-lived P. falciparum (Carter and Mendis 2002).

A powerful piece of circumstantial evidence supports the existence of relatively longstanding P. vivax infection in pre-agricultural Africa: Duffy antigen negativity. The Duffy antigen is a chemokine expressed on RBC membranes, and also happens to be an essential receptor for P. vivax merozoite entry into RBCs (Carter and Mendis 2002). Homozygotes for Duffy negativity are thus completely immune to P. vivax, and moreover appear to suffer no ill health-effects. In native populations, Duffy negativity prevalence is almost 100% in most west and central Africa (Culleton and Carter 2012), likely the ancestral seat of malaria and the areas of the most intense malaria transmission today, while Duffy negativity is highly prevalent throughout the rest of the continent. Since homozygosity is required for significant benefit, one may that expect tens of thousands of years are necessary for Duffy negativity to become fixed in a population (Carter and Mendis 2002), and Hamblin and Rienzo (2000) estimated a selective sweep may have occurred 33,000 years ago (95% CI 65,000–97,200 years ago). Thanks to Duffy negativity, P. vivax was likely driven nearly to extinction in Africa, but escaped into Asia and the larger world (Liu et al. 2014), perhaps around 10,000 years ago, where populations have had insufficient time to evolve Duffy negativity (Culleton and Carter 2012).

About 10,000 years ago, African proto-agriculture led to more sedentary, larger human settlements that could sustain more virulent, short-lived infections (Carter and Mendis 2002). It was around this time that P. falciparum in gorillas may have crossed over into humans, according to a recent hypothesis (Loy et al. 2017). Even if the gorilla hypothesis is false, it is clear that P. falciparum did not become a significant human disease until between 5000 and 10,000 years ago, and that its rise was related to that of agriculture (Carter and Mendis 2002; Webb 2014). Malaria, both P. vivax and P. falciparum, likely spread through most of the inhabited world during early historical times (i.e. before the common era), although P. vivax mainly affected the more northerly regions, given its dormant phase and better cold tolerance versus P. falciparum (Packard 2007); malaria was rapidly introduced to the New World following its discovery by Europeans.

In the nineteenth century, malaria reached its global zenith, with most of the globe’s population at risk (Carter and Mendis 2002), but then declined into extinction in most of Europe and the Americas by the mid-twentieth century, its retreat primarily caused by agricultural modernization and changing living conditions that discouraged transmission, and aided by later eradication programs (Packard 2007; Carter and Mendis 2002). This, however, was not the experience of tropical Africa, and from here out we will restrict our attention to this continent.

3.3 The colonial era, Africa, and modern malariology

In the late 1800s, spurred largely by the discoveries of Koch and Pasteur, the search was on for bacterial causes of many diseases, and in 1880, Charles Laveran, an obscure French army officer stationed in Algeria (a French colony at the time, having been subdued in a bloody war of conquest spanning 1830–1847, and that killed as much as a third of the native population (Kiernan 2007)), observed a variety of strange writhing forms within the erythrocytes of malaria victims, which he would come to identify as a protozoan parasite that he named Oscillaria malariae. It was the first protozoan discovered to infect man, and Laveran would receive the Nobel Prize in 1907 for this discovery (Cox 2010).

It fell chiefly to Ronald Ross, a British physician, to elucidate the vector by which the malarious protozoan was transmitted, the female anopheline mosquito, which he showed in birds in 1897, and in humans in Freetown, Sierra Leone, in 1899 (Cox 2010). Note that while to modern ears, the idea of a mosquito transmitting a disease is entirely natural, at that time it was a truly novel notion (Cox 2010). Sierra Leone had been established as a British colony in 1787, and due to the high malaria burden came to be known as the “White Man’s Grave” (Bockarie et al. 1999), and indeed, malaria has been credited by some historians as protecting the interior of Africa from the depredations of European colonialism during the slave era (Webb 2014). The discovery of infected anopheline vectors A. gambiae and A. funestus, along with their breeding sites in myriad small pools, by Ross and his colleagues during their 1899 expedition led to vector control measures including bednets, window screens, and larval control via oiling of pools (Bockarie et al. 1999; Webb 2014).

Ross also recommended segregating European and African populations to protect the Europeans (Bockarie et al. 1999); this too would be a feature, although varying in degree by time and place, of colonial malaria control efforts (Webb 2014).

In Sierra Leone and elsewhere, subsequent efforts included draining or oiling pools, and removing household receptacles that could support breeding. Other anti-larval efforts, of which Ross was a champion, included treating pools with a highly toxic copper-based compound known as Paris Green, stocking with larvivorous fish, and, by World War II, treating with oils containing the pesticide DDT (dichloro-diphenyl-trichloroethane) (Bockarie et al. 1999; Webb 2014). While sometimes effective, anti-larval measures required ongoing action, were labor-intensive, and dependent upon funding. A common pattern was concentrating malaria control efforts in urban areas, and in commercial areas where European interests desired a healthy indigenous workforce, but with European health as a priority. There was also legitimate concern that measures decreasing malaria prevalence among native populations could reduce immunity, rendering them vulnerable to epidemic malaria (Webb 2014).

DDT, first used against malaria by the US Army in World War II, has a long-lasting residual effect, such that a dwelling need be sprayed only infrequently to have a toxic effect on mosquitoes. Thus, the 1940s and 50s ushered in the pesticide era, with indoor residual spraying (IRS) increasingly used by national control programmes (Nájera et al. 2011). Macdonald’s mathematical model (Macdonald 1957) (discussed in Sect. 4.2) provided a powerful theoretical basis for increasing adult mortality as the linchpin of control (Nájera et al. 2011), and furthermore, the strategy of spraying was viewed as general and inexpensive, compared to expensive quinine treatment or labor- and capital-intensive environmental engineering and larvaciding (Webb 2014). Against this background, the WHO launched its Global Malaria Eradication Programme (GMEP) in 1955, based on spraying with DDT and related compounds supplemented by mass drug treatment, and efforts were geared toward eradication over control, with malaria control viewed with contempt by the program’s architects (Nájera et al. 2011). The GMEP coordinated with national control programmes, and launched a number of pilot projects, the most famous and well-done being the Garki Project, which motivated the Garki mathematical model, by Dietz et al. (1974) (Sect. 4.4).

A deeply unfortunate effect of the GMEP was the undermining of the specialized field of malariology (why study malaria if all one need do is spray DDT, regardless of the particular vector, geography, socially, biology, etc.?) (Nájera et al. 2011; Webb 2014), as well as the (temporary) abandonment of many control measures other than IRS (Nájera et al. 2011). Despite very meaningful successes, including the eradication of malaria from many countries (especially outside of Africa), the GMEP was beset by setbacks, and in 1969 it was acknowledged that eradication was not a realistic short-term goal in many regions, marking the end of the programme (Nájera et al. 2011).

Beginning in the 1960s, the synthetic anti-malarial chloroquine became widely and inexpensively available on the African continent, with dramatically positive health consequences, being chiefly responsible for marked reductions in child and malaria-specific mortality through the 1960s and 70s (Carter and Mendis 2002; Webb 2014). However, malaria began a resurgence throughout Africa beginning in the late 1970s and 1980s, largely attributable to the evolution and spread of chloroquine-resistant P. falciparum (Carter and Mendis 2002). However, financial- and debt-crises borne by the recently independent African states, reductions in public health spending, widespread and large-scale political violence and chaos, the HIV/AIDS epidemic, and the widespread expansion of rural agriculture into west and central African rainforests, where deforestation created new habitat for anopheline vectors, all played roles (Webb 2014). It was also revealed in this era that prior concerns that “protecting” populations in endemic zones where elimination was infeasible could dangerously undermine immunity were well-founded, as deadly epidemics swept through many such regions, most notably in the highlands of Madagascar in 1986 (Carter and Mendis 2002; Webb 2014).

In the face of devastating infectious disease across the Global South, the WHO and several other organizations founded the Roll Back Malaria Partnership in 1998, while the WHO’s “Global Fund to Fight AIDS, Tuberculosis and Malaria” was established in 2002, and in 2007, the Bill and Melinda Gates Foundation announced a campaign to eradicate malaria (Webb 2014). New tools became available, mainly insecticide-treated bednets (ITNs), and the burden of chloroquine resistance was relieved with newer artemisinin-based combination therapy. These campaigns have enjoyed apparent success, with a 57% decrease in African malaria mortality (per 10,000) from 2000 to 2015 (Gething et al. 2016); Bhatt et al. (2015) estimated that 68% of avoided malaria cases in Africa (from 2000 to 2015) could be attributed to ITNs.

Whether these gains will be maintained has yet to be seen. Malaria eradication and control programs have historically seen their greatest success in the first few years (Webb 2014), not all countries have experienced similar improvements under similar control programs (Snow et al. 2015), malaria incidence has recently increased locally in some areas, e.g. coastal Kenya (Snow et al. 2015), Plasmodium artemisin resistance has emerged in southeast Asia (Webb 2014), and perhaps even more worrisome, widespread pyrethroid resistance (the insecticide in ITNs) is evolving in vectors across Africa (Hemingway et al. 2016), although the impact of these developments in resistance is yet to be proven. Furthermore, a general dynamic of increased control of a childhood illness is an early drop in disease transmission, and a consequent shift in disease burden from younger to older ages that generates a rebound increase in incidence after a few years (Griffin et al. 2016). In the context of malaria, an in-depth modeling study by Griffin et al. (2016) suggested that merely sustaining current control efforts, even absent new vector or parasite resistance, will lead to increases in malaria incidence and mortality by 2020. And finally, climate change continues in its insidious trajectory, with uncertain consequences.

4 Early mathematical models of malaria

4.1 Sir Ronald Ross, a pioneer of quantitative epidemiology

Ross proposed several mathematical models for malaria transmission in the early 1900s that were analyzed and modified by others, including Alfred J. Lotka (Smith et al. 2012). These were extended by Macdonald and other investigators in the 1950s, and this work, which focused heavily on ${\mathcal R}_0$ and mosquito eradication for malaria control, would prove to be very influential in guiding the ultimately failed GMEP (1955–1969). Smith and colleagues, who expertly reviewed the early history of Ross and Macdonald (Smith et al. 2012), have argued that there is no single or canonical “Ross–Macdonald” model, and that it is more instructive to understand this as a family of models characterized by a set of broadly shared assumptions and key entomological and epidemiologic parameters, whose estimation was historically motivated by quantitative models. Note, however, there is a clear “Macdonald” model, as presented in Macdonald (1957).

Ross’s original 1908 model is of purely historical interest, so we will skip to the 1911 model (which was solved and extensively analyzed by Lotka) given as (Smith et al. 2012)

$$\begin{aligned} \frac{dX}{dt}= & {} m a b z (H - X) - r X, \end{aligned}$$

(1)

$$\begin{aligned} \frac{dZ}{dt}= & {} a c x (M - Z) - g Z, \end{aligned}$$

(2)

where H is the total human population with X the infected component, M and Z are similarly the total and infected mosquito populations, m is M / H (mosquitoes/man), a is the mosquito biting rate (bites/mosquito/day), z is the infectious mosquito fraction (Z / M), r is the human recovery rate (day$^{-1}$), b is the probability of human infection after an infectious bite (omitted and implicitly 1 in Ross’s original model, we include it for clarity), c is the probability of a human infecting a mosquito upon biting, $x = X/H$ is the parasite rate, and g is the mosquito death rate (day$^{-1}$).

Sharpe and Lotka extended this model, as reviewed by Smith et al. (2012), to include latency between inoculation and infectivity in both man and mosquito, but because their model failed to consider mosquito mortality during latency, biological conclusions were flawed.

4.2 Ross–Macdonald

In a series of works in the early 1950s, Macdonald formulated a highly influential model, based on Ross’s basic model. The major mathematical innovations introduced by Macdonald over Ross were accounting for the delay to infectiousness in mosquitoes, and superinfection, where multiple malarial strains can coinfect a host and are independently cleared, thus altering the recovery rate from the infected to recovered/susceptible category. Unfortunately, this concept was incorrectly translated into mathematics by Macdonald, apparently due to a miscommunication (Smith et al. 2012). The correct form was described by Walton in 1947 (per Dietz et al. (1974)), and was incorporated into the influential “Garki” model devised by Dietz et al. in 1974 (Dietz et al. 1974). We have, using x(t) as the proportion of infected humans (the “parasite rate”), Macdonald’s model (Macdonald 1957; MacDonald et al. 1968)

$$\begin{aligned} \frac{dx}{dt} = h (1 - x) - \rho (r,h) x, \end{aligned}$$

(3)

where h is the inoculation rate and r is the first-order rate of recovery from each infecting malarial strain, each of which is assumed to be cleared independently; the overall rate of recovery, $\rho (r,h)$, is a function of the inoculation rate and strain-specific recovery rate. Now, inoculation is given as

$$\begin{aligned} h = \frac{m a^2 b c p^n x}{a x - \ln (p)} = \frac{m a^2 b c x}{a x + g} \exp (-n g), \end{aligned}$$

(4)

where n is the duration of the sporogonic cycle, p is the daily probability of survival, implying $p = \exp (-\,g)$ and that $\exp (-\,n g)$ is the fraction of mosquitoes surviving from the time of exposure to infectivity; other parameters are as in the Ross model (c was assumed to be unity by Macdonald (1957), but we have included it for generality). In his 1957 book (Macdonald 1957), Macdonald derives this expression as

$$\begin{aligned} h = m a b c s, \end{aligned}$$

(5)

where s is the sporozoite rate (i.e. the fraction of mosquitoes with sporozoites in their salivary glands) which in turn is derived as follows. We have, from the exponential distribution, that the expected (mean) lifetime of any mosquito is

$$\begin{aligned} \frac{1}{g} = \frac{1}{-\ln (p)}. \end{aligned}$$

(6)

We then have that the total mosquito-days spent in a potentially infectious state, i.e. they have survived at least n days (duration of the sporogonic cycle), is

$$\begin{aligned} \frac{p^n}{-\ln (p)}. \end{aligned}$$

(7)

To determine the sporozoite rate, we need know what fraction of such potentially infectious mosquitoes actually are infectious. We have that the average number of infectious feeds in a day is ax, and so the probability of taking no infectious feeds in a day is $\exp (-ax)$, and the chances of both surviving and taking no infectious feeds is $p \exp (-ax)$. It follows (again from the exponential distribution) that the expected life taking no infectious feeds is

$$\begin{aligned} \frac{1}{-\ln \left( p \exp (-a x)\right) } = \frac{1}{ax - \ln (p)}, \end{aligned}$$

(8)

and the total mosquito-days spent in the potentially but non-infectious state is

$$\begin{aligned} \frac{p^n}{ax - \ln (p)}. \end{aligned}$$

(9)

We finally arrive at the sporozoite rate (i.e. fraction of mosquitoes in an infectious state) as

$$\begin{aligned} s = \frac{(7) - (9)}{(6)} = \frac{p^n a x}{ax - \ln (p)}. \end{aligned}$$

(10)

This expression can also be derived by applying a quasi-steady-state assumption to a delay-differential version of Ross’s model,

$$\begin{aligned} \frac{dx(t)}{dt}= & {} m a b z(t) (1 - x(t)) - \rho (r,h) x(t), \end{aligned}$$

(11)

$$\begin{aligned} \frac{dz(t)}{dt}= & {} a c x(t - n) (1 - z(t - n))\exp (-n g) - g z(t). \end{aligned}$$

(12)

That is, setting dz / dt = 0 (assuming $x(t - n) = x(t)$) and solving for z. Moving on, the recovery rate, $\rho (r,h)$, takes the form

$$\begin{aligned} \rho (r,h) = \left\{ \begin{array}{lr} r - h, &{} h < r \\ 0, &{} h \ge r \end{array} \right. , \end{aligned}$$

(13)

but this implies no recovery ever occurs when inoculation exceeds the strain-specific recovery rate (clearly an error). The correct form, given by Dietz et al. (1974), is

$$\begin{aligned} \rho (r,h) = \frac{h}{\exp \left( \frac{h}{r}\right) - 1}. \end{aligned}$$

(14)

Finally, from Macdonald’s model, we can derive the following expression for ${\mathcal R}_0$:

$$\begin{aligned} {\mathcal R}_0 = \frac{m a^2 b c p^n}{-r \ln (p)} = \frac{m a^2 b c \exp (-n g)}{r g}. \end{aligned}$$

(15)

The key conclusion from this expression is that daily mosquito survival, p, appears in both the numerator and the denominator, such that decreasing it lowers ${\mathcal R}_0$ logarithmically (Macdonald’s ${\mathcal R}_0$, as a function of several different parameters of the first equality in Eq. (15), is depicted in Fig. 5). This suggests targeting the adult mosquito vector as the most efficacious strategy for malaria control, and indeed, this was the basic approach of the GMEP, which relied principally upon indoor residual spraying with DDT (and other pesticides) to achieve this end, as discussed already in Sect. 3.3.

While ${\mathcal R}_0$ in Eq. (15) suggests targeting mosquito daily survival, p, over mosquito density, m, it is obvious from the basic ecology that these are not independent parameters. Nor, indeed, is a, the biting rate, since biting provides blood needed to nourish the mosquito’s eggs. Both Ross and Macdonald were pioneering thinkers, but it seems clear that a more robust model framework that more fully accounts for the vector lifecycle is necessary for us to be confident in any conclusions. We shall explore some of these issues in detail in Sect. 6.

Fig. 5

Change in Macdonald’s ${\mathcal R}_0$ as a function of each major parameter (except b and c, which have a straightforward linear effect). Daily mosquito survival, p, is the most sensitive parameter, and ${\mathcal R}_0$ is also given as a function of expected mosquito life, $1/\ln (p)$

4.3 Vectorial capacity

In 1964, Garrett-Jones (1964) proposed an alternative metric, contra ${\mathcal R}_0$, for assessing and motivating vector control, namely the vectorial capacity (VC). It was defined qualitatively, for a vector population, as “the average number of inoculations, ..., originating from one case of malaria in unit time [typically in days], that the [vector] population would distribute to the human host if all female adult mosquitoes biting the human host became infected” (Garrett-Jones 1964). In other words, it is the number of new malaria cases (i.e. infectious bites, assuming all such bites result in infection) originating from a single case in a single day. Garrett-Jones (1964) formally defined it as the product of (1) the man-biting rate (total bites/person/day), which is the total number of mosquitoes infected from a single case in a single day, (2) the expectation of infective life, and (3) the man-biting habit, the number of bites on the human host per day per individual mosquito. Using the Ross–Macdonald parameters, the vectorial capacity is given by Garrett-Jones and Shidrawi (1969)

$$\begin{aligned} \text {VC} = \frac{m a^2 c p^n}{-\ln (p)} = \frac{m a^2 c \exp (-n g)}{g}. \end{aligned}$$

(16)

Note that the original form implicitly assumed c = 1, but we have included it for generality. The vectorial capacity is very similar to the concept of ${\mathcal R}_0$, but is simply considering the number of new cases (assuming 100% infectiousness of bites) that result in the first unit of time from the original case, rather than over the lifetime of the first case. Where ${\mathcal R}_0$ has units of cases, VC has units cases/day, and the two terms relate (again, under Macdonald’s model) as

$$\begin{aligned} {\mathcal R}_0 = \text {VC} \times \underbrace{b}_{\text {probability that man is infected by a bite}} \times \underbrace{\frac{1}{r}}_{\text {Expected time that first case is infectious}}. \end{aligned}$$

(17)

The concept of vectorial capacity was used in a number of subsequent quantitative studies, such as Garrett-Jones and Shidrawi (1969), Dietz et al. (1974), Molineaux et al. (1978), and it is also noteworthy that VC is the component of ${\mathcal R}_0$ that is most directly affected by weather parameters (Craig et al. 1999).

4.4 Developments post-Ross–Macdonald

The next major mathematical modeling contribution to malaria transmission dynamics was by Dietz et al. (1974) and entailed the inclusion, into the basic Ross–Macdonald framework, of a kind of slowly-acquired immunity that results in a non-infectious parasitemia following inoculation that is cleared relatively rapidly. Non-immune hosts are assumed to manifest infectious clinical disease that transitions to a non-infectious parasitemia that is cleared slowly. The model was fit to data for two villages in the Garki district of Nigeria, where data on the parasite and sporozoite rate had been collected by age. A follow-up work by Molineaux et al. (1978) in 1978 compared the model against several other datasets, where vectorial capacity was estimated from entomologic parameters and observed host-biting rates.

Many of the major modeling contributions following this work concern the proper or realistic modeling of immunity, especially the distinction between anti-disease (resistance against the harmful clinical manifestations of parasite infection, such as fever, anemia, etc.) and anti-parasite immunity (resistance against the actual Plasmodium infection), and how these are induced with infection and lost with time. Since climate, and not immunity, per se, is our primary focus, we defer a brief discussion of this model genealogy to Sect. 7.1. Moreover, most climate-focused models have only included fairly rudimentary descriptions of immunity, if it is addressed at all (but see (Yamana et al. 2013, 2017) for exceptions), and the hybridization of these two modeling traditions is a major future challenge.

5 Quantifying the relationships between weather and the parasite and vector lifecycles

Understanding the quantitative relationships between weather, primarily temperature and rainfall (and to a lesser degree, relative humidity), and the malaria parasite and vector lifecycles is critical to a realistic and meaningful assessment of the impact of current and projected climate change on malaria transmission dynamics. We review some of the most widely used quantitative relationships and the data they are based here, but it should be understood that these data are drawn from a variety of Plasmodium and Anopheles species, and the widely used formula of Moshkovsky for sporogonic and gonotrophic durations, for example, is based on data from the 1930s obtained in the European vector A. maculipennis (Detinova 1962). Moreover, there is some suggestion that thermal sensitivities may vary between laboratory and field strains of P. falciparum, which may be more adapted to local conditions, although Lyons et al. (2012) observed similar thermal tolerances between laboratory and wild strains A. arabiensis and A. funestus. While still poorly understood in general, short- and long-term adaptations in both vector and parasite to local temperature ranges and shifts could limit the ultimate generalizability of model inferences made assuming thermal response functions uniform throughout space and time (Sternberg and Thomas 2014).

Furthermore, while only a single thermal response function for a given process is typically considered in models, major vectors may differ importantly in how they respond to both mean and fluctuating temperatures, with Lyons et al. (2013) observing survival and development in the three major African vectors, A. gambiae, A. arabiensis, and A. funestus to vary both in response to mean temperature and temperature fluctuations, with A. funestus in particular much more sensitive to temperature fluctuations than A. arabiensis. Moreover, most models use thermal response functions drawn from multiple species, and thus how explicit consideration of the varying responses to weather between relevant vectors may affect model conclusions remains an open question.

5.1 Parasite

The sporogonic cycle of Plasmodia, i.e. infection and sexual reproduction in the mosquito midgut to ultimately yield infectious saliva sporozoites, is very clearly directly influenced by climate, with warmer temperatures (at least to a point), leading to more rapid parasite development (decreasing n, under the Ross–Macdonald framework), and this has been the focus of most mathematical works. We also note, however, that temperature may affect infectivity to both mosquito (c, per Ross–Macdonald) and man (b). Temperatures above about 30 $^{\circ }\hbox {C}$ may decrease P. falciparum survival in the mosquito midgut, and hence decrease c (Eling et al. 2001; Okech et al. 2004a), and Paaijmans et al. (2012) (in a rodent model) found increasing temperatures to decrease the prevalence of sporozoites in infected mosquitoes, and hence decrease b. However, these factors are less frequently accounted for in models, and we restrict further attention to the sporogonic cycle and n.

5.1.1 Sporogonic cycle and temperature

The classical formula of Moshkovsky It has long been recognized that the duration of the sporogonic (extrinsic) cycle in mosquito, denoted by n, is hyperbolically related to temperature. That is, given a constant D, measured in degree-days (the “sum of heat,” as elaborated below), a minimum temperature, $T_{min}$, and mean ambient temperature, $T > T_{min}$ (in $^{\circ }\hbox {C}$), we have

$$\begin{aligned} n = \frac{D}{T - T_{min}}. \end{aligned}$$

(18)

This relation is based on a now obscure 1935 work in Russian by Nikolaev (1935), who, as related by Detinova (see pp. 122–150 of Detinova (1962)), gathered data on the duration of sporogony in the A. maculipennis mosquito, the historic vector of malaria in Europe and the Middle East (Djadid et al. 2007), and one that has recently been found to be expanding its range northeasterly into Eurasia, likely as a consequence of global warming (Novikov and Vaulin 2014). Nikolaev (1935) determined that P. falciparum sporogony ceases below 18 $^{\circ }\hbox {C}$ (16 $^{\circ }\hbox {C}$ for P. vivax), and his data was used by later authors (Moshkovsky in particular) to estimate the natural length of the sporogonic cycle.

Briefly, the Moshkovsky method (as given in Detinova (1962)) assumes the “sum of heat” hypothesis that a fixed amount of heat—which can be measured in degree-days calculated beyond the minimum temperature required for sporogony to progress at all—is required to complete the cycle. This “sum of effective temperatures” was estimated at 105 $^{\circ }\hbox {C}$-days for P. vivax, 111 $^{\circ }\hbox {C}$-days for P. falciparum, and 144 $^{\circ }\hbox {C}$-days for P. malariae. Given daily mean diurnal temperature, the time to completion of sporogony can be straightforwardly calculated for any day of the year. The formula of Moshkovsky has been used in many more modern works (Molineaux et al. 1978; Craig et al. 1999; Hoshen and Morse 2004; Parham and Michael 2010). It should, however, be emphasized that, even in the early days, it was observed that temperatures above 32 $^{\circ }\hbox {C}$ were lethal to developing parasites, blocking sporogony (Macdonald 1957). Sporogonic duration as a function of constant temperature is given in Fig. 6.

Recent developments More recent studies are consistent with lower temperatures prolonging sporogony, but several laboratory works suggested (constant) temperatures of only 30 $^{\circ }\hbox {C}$ could significantly impede early P. falciparum sporogony in A. stephensi (a vector commonly found in the Indian subcontinent) (Noden et al. 1995; Eling et al. 2001). This finding’s generality was challenged, however, by work by Okech and colleagues (Okech et al. 2004a, b) who used wild P. falciparum strains in western Kenya. They found that naturally fluctuating field temperatures, up to 33 $^{\circ }\hbox {C}$ maximum, did not interrupt parasite development in A. gambiae (Okech et al. 2004b), and under laboratory conditions constant exposure to 32 $^{\circ }\hbox {C}$ decreased, but did not eliminate, early P. falciparum survival in mosquito midgut and infectivity (Okech et al. 2004a). These results indicate that wild P. falciparum typically exposed to hotter conditions is more thermally adapted than lab-raised strains, and thus it seems likely that the Moshkovsky formula amended by sporogonic arrest around 32–34 $^{\circ }\hbox {C}$ is a reasonable description for wild P. falciparum, although mosquito infectivity may decrease at 30–32 $^{\circ }\hbox {C}$ (Okech et al. 2004a). Alternatively, several authors (Paaijmans et al. 2009; Mordecai et al. 2013) have more recently used a general unimodal functional form given by Briere et al. (1999),

$$\begin{aligned} r(T) = c T (T - T_0) (T_m - T)^{\frac{1}{2}}, \end{aligned}$$

(19)

with Mordecai et al. (2013), for example, using $\hbox {c} = 0.000111$, $T_m = 34.4~{}^{\circ }\hbox {C}$, and $T_0 = 14.7~{}^{\circ }\hbox {C}$. Such a formulation gives developmental rates qualitatively similar to those under amended versions of Moshovsky’s formula, as shown in Fig. 6.

Fig. 6

Duration (left panel) and rate (right panel) of the sporogonic cycle as a function of constant temperature, per Moshkovsky’s formula for P. falciparum (D = 111, $T_{min}$ = 16) and P. vivax (D = 105, $T_{min}$ = 14.5) (Detinova 1962), and assuming sporogony ceases above about 33 $^{\circ }\hbox {C}$. Curves calculated from Briere’s formula, using parameter values from Mordecai et al. (2013), are also given for comparison

5.2 Vector

5.2.1 Gonotrophic cycle, oviposition, and biting rates

Adult female Anopheles mosquitoes take blood meals from human hosts almost exclusively to provide energy and nutrients for their eggs, and the gonotrophic cycle (defined as the time between blood meals) is classically divided into the following three stages (Detinova 1962):

• 1. Search for host and attack.

• 2. Digestion of blood meal and egg maturation. This stage is highly temperature-dependent.

• 3. Search for body of water and oviposition. This stage is dependent upon the availability of suitable standing waters, itself a function of recent rainfall (except for those vectors that breed in more permanent bodies of water).

The terminology around gonotrophy is sometimes abused, with either stage II (blood meal digestion) or stages II and III (blood meal digestion and oviposition) together sometimes conflated with the gonotrophic cycle as a whole. Gonotrophic cycle duration is important because, firstly, it partially determines mosquito population levels. Second, it has been traditionally assumed that mosquitoes take only a single blood meal between ovipositions (Scott and Takken 2012), implying that the duration of the gonotrophic cycle is identical to the biting interval. There is a biochemical basis for this assumption, with mosquito host-seeking behavior inhibited following a blood meal via two mechanisms that act in sequence (Takken et al. 2001). First, gastric distention in the period immediately following feeding activates inhibitory mechanoreceptors. Second, blood digestion and adequate nutritional status seem to stimulate neuropeptides that decrease the sensitivity of olfactory receptors in the antennae (Takken et al. 2001). Evolutionarily, avoiding unnecessary blood meals has obvious potential benefits, as biting requires energy to seek out hosts who, as we might expect, are then rather hostile to blood-sucking insects (Klowden and Briegel 1994), although as a nocturnal feeder biting somnolent hosts, A. gambiae may have faced less such selective pressure (Klowden and Briegel 1994).

Multiple blood feeding The assumption of one blood meal per gonotrophic cycle is clearly violated for at least a subset of A. gambiae mosquitoes, namely newly-emerged adult females with insufficient nutritional reserves to fuel a gonotrophic cycle from a single blood meal (Scott and Takken 2012). These mosquitoes, whose ovaries do not progress beyond Christophers’ stage II (an early stage in ovary development) after their first blood meal are termed “pre-gravid,” and typically take a second blood meal within 24 h to stimulate egg development (Scott and Takken 2012). Unsurprisingly, it is smaller mosquitoes that appear to be more nutritionally deficient, requiring two blood meals to complete their first gonotrophic cycle, with the first blood meal apparently devoted mainly to shoring up nutritional reserves, while larger individuals can successfully complete oogenesis with a single blood meal (Takken et al. 1998). After this early transient phase, surviving small and large adult mosquitoes may behave similarly, but with smaller mosquitoes now more likely to be infected by Plasmodium (Takken et al. 1998). Emerging adult mosquito size is dependent on larval conditions, with higher densities, higher temperatures (which increase the rate of larval development), and poor nutrition favoring smaller adult sizes (Smith et al. 2012). It should then be noted that temperature could have a subtle second-order effect on malaria transmission by increasing the proportion of smaller adult female mosquitoes more likely to engage in multiple biting behavior. Further, most “extra” bites would occur early in life before the completion of sporogony.

Whether multiple blood feeding is common among anophelines outside of this context is unclear. Klowden and Briegel (1994) observed laboratory A. gambiae to strongly seek human hosts every 24 h following a blood meal, corresponding to regular nocturnal feeding behavior and suggesting no inhibition of host-seeking. In direct contrast, however, in laboratory-raised A. gambiae of uniform size (thus avoiding contamination by smaller pre-gravid mosquitoes), Takken et al. (2001) observed blood meals to strongly inhibit host seeking up to 72 h following a blood meal (at 27 $^{\circ }\hbox {C}$), with this interval corresponding to egg maturation or oviposition in most individuals. Despite its potential importance to malaria epidemiology, the pre-gravid anopheline is typically neglected in mathematical models, although Depinay et al. (2004) took the probability of taking multiple bloodmeals to be a function of both mosquito weight and parity.

Quantifying the gonotrophic cycle Leaving the complications of possible multiple blood feedings aside, let us now review quantifications of the gonotrophic cycle duration, or (typically) equivalently, biting rate. The classic works of Macdonald (Macdonald 1956b, 1957) and other derived works (e.g., (Molineaux et al. 1978)) have generally assumed a constant biting interval of about two days for A. gambiae (Molineaux et al. (1978) chose two days on the basis of the distribution of abdominal stages in sprayed mosquitoes collected in the Garki Project), although some older field studies suggested three days as more likely (Takken et al. 2001). However, it has been well-understood since at least the 1950s that stage II varies greatly in inverse proportion to temperature (Macdonald 1957).

Duration of stage II Detinova (1962), in his influential work, presented data on the duration of stage II, which also followed the basic form of Moshkovsky’s degree-day formula, but the constants varied with relative humidity (RH), such that this stage is quicker under more humid conditions; curves and formula constants are given in Fig. 7. More recent studies include (Ra et al. 2005; Afrane et al. 2005; Lardeux et al. 2008; Mala et al. 2014). Probably most notable is a 2008 laboratory study (Lardeux et al. 2008) on A. pseudopunctipennis, a major South American vector, who studied the duration of stage II/III (i.e. blood meal to oviposition) under a range of temperatures, and used the following general function, adopted from Lactin et al. (1995) (based in turn on Logan et al. (1976)) for development rate, r(T), assumed to be the reciprocal of the mean time to oviposition (i.e., it is assumed that gonotrophic duration is exponentially distributed):

$$\begin{aligned} r(T) = \exp (\rho T) - \exp \left( \rho T_m - \frac{(T_m - T)}{\varDelta }\right) + \lambda , \end{aligned}$$

(20)

where $T_m$ is the “thermal maximum,” or lethal temperature, $\varDelta$ is the temperature range over which the development rate begins to fall from a maximum to zero (“temperature boundary layer”), $\rho$ determines the increase in development rate with temperature and, per Logan et al. (1976), can be “interpreted as a composite value for critical enzyme-catalysed biochemical reactions,” and lastly, $\lambda$ gives r(T) when T = $T_m$. Development time reached its nadir at 31 $^{\circ }\hbox {C}$, and there was no oviposition at 37 $^{\circ }\hbox {C}$, as all mosquitoes suffered mortality before they could oviposit.

Now, a subtlety here is that a single development rate does not fully describe the data, and because mosquitoes lay their eggs only at night, the actual time to oviposition clusters at intervals of 24 h. For example, the overall mean time to oviposition at 35 $^{\circ }\hbox {C}$ was 2.3 days, but roughly 70% of mosquitoes oviposited on night two, and 30% on night three. The lower the temperature, the longer the mean time, and oviposition events are spread over more nights; this is shown in Fig. 7, where simulated cohorts of ovipositing mosquitoes are generated using parameters fit to Eq. (20) as reported in Lardeux et al. (2008). It is obvious from this data that time to oviposition is not exponentially distributed, and the epidemiological implications could be addressed with some future model.

Duration of stages I and III Work by Detinova (1962) from 1953 on A. maculipennis in the USSR, in which the contraction of ovarioles (which begins following oviposition) was determined in dissected mosquitoes, suggested that a majority of mosquitoes take a blood meal within 8 h of oviposition, and over two-thirds of mosquitoes have taken a blood meal within 24 h of oviposition. With respect to the third stage, this author also determined that oviposition usually occurs within 24 h of egg maturation. Thus, it is reasonable to assume that the third and first phase take, together, around 24 h (but occasionally up to 48 h), while the second phase is highly temperature-dependent. Note, however, that low availability of oviposition sites (related, say, to low rainfall) could significantly prolong search time, as studied in a simple model by Gu et al. (2006).

Overall duration The overall duration of the gonotrophic cycle can be reasonably estimated as 24 h for stage I and III combined, plus a temperature dependent term for stage II, either Moshkovsky’s formula, or a relation derived from other data, such as that of Lardeux et al. (2008). These two options (and using several relative humidities for Moshkovsky’s formula) are compared in the right panel of Fig. 7.

Fig. 7

The left panel shows simulated proportions of mosquitoes that oviposit, as a function of time, at different ambient temperatures based on the experiments of Lardeux et al. (2008). Oviposition is spread over several nights, with the number of nights increasing as mean time to oviposition increases. The blue histograms show simulated time to oviposition if all events occur on a single night (curve is scaled down for clarity). The right panel compares the duration of the gonotrophic cycle as measured by Lardeux et al. (2008), and as calculated by Moshkovsky’s formula plus an additional day for stages I and III, and using, per Detinova (1962), $D = 65.4$, $T_{min} = 4.5$; $D = 36.5$, $T_{min} = 9.9$; and $D = 37.1$, $T_{min} = 7.7$, for relative humidities (RH) of 30–40%, 60–70%, and 90–100%, respectively

Time to first bloodmeal We finally note that the there is a small delay, on the order of 1–3 days, between emergence from pupal stage to the first bloodmeal (“pre-bloodmeal period”) and the start of the gonotrophic cycle proper (Paaijmans et al. 2013a). Furthermore, Paaijmans et al. (2013a) found this pre-bloodmeal period to be temperature-dependent in An stephensi, being about three days at 18 $^{\circ }\hbox {C}$ but just one day at either 26 or 32 $^{\circ }\hbox {C}$.

Thus, we have the pre-bloodmeal period as an additional adult stage that is typically disregarded in models, but acts to delay time to first infection and infectivity in a temperature-dependent manner. This phenonemon may slightly narrow the temperature range over which malaria is effectively transmitted (Paaijmans et al. 2013a), and moreover, will shift infectivity to older mosquito age classes, which could in turn interact with age-dependent temperature-mediated mortality (see Sect. 5.2.2).

5.2.2 Daily mosquito survival

Age, Plasmodium infection, and temperature all influence mosquito survival. Classically, it has been assumed that random mortality events such as predation, etc. account for most mosquitoes deaths, such that senescence can be disregarded and death modeled as a Poisson process, i.e. survival times are exponentially distributed (this is also the implicit assumption of those mathematical models that employ first-order death kinetics in differential equations), and daily survival probability is typically assumed to be around 0.90–0.95 (Macdonald 1956a; MacDonald et al. 1968). Nevertheless, mosquito death hazard clearly increases with age (Clements and Paterson 1981; Okech et al. 2003; Dawes et al. 2009; Christiansen-Jucht et al. 2014), with this most convincingly demonstrated in laboratory studies (e.g., Bayoh 2001; Christiansen-Jucht et al. 2014), although Clements and Paterson (1981) found evidence of age-dependent mortality in wild populations as early as 1981, and recent analysis by Ryan et al. (2015a) also suggested senescence occurs in wild mosquitos, despite high extrinsic mortality rates. Survival is often described by the Gompertzian distribution (Clements and Paterson 1981; Bayoh 2001; Christiansen-Jucht et al. 2014; Ryan et al. 2015a), given as

$$\begin{aligned} S(t; \lambda , \theta ) = \exp \left( \frac{\lambda }{\theta } \left( 1 - \exp (\theta t)\right) \right) , \end{aligned}$$

(21)

with the model essentially positing exponentially increasing mortality with age, although other mathematical descriptions (e.g., double exponential model, quadratic model) can also describe age-dependent mortality increases (Clements and Paterson 1981), and one especially useful description is the gamma distribution, which may be straightforwardly implemented in an ordinary differential equations (ODE) setting, as discussed further in Sect. 6.3.5.

Laboratory mortality is straightforward to measure; wild mosquito survival can be estimated by means of mark-release-recapture (MRR) experiments. Survival time is, in this case, often fit to an exponential model, as in, e.g. (Midega et al. 2007; Olayemi and Ande 2008), but see Ryan et al. (2015a) for a recent example of applying the Gompertz distribution to wild populations.

Temperature and humidity dependence Temperature strongly affects mosquito survival, with (age-dependent) death rates increasing above about 23 $^{\circ }\hbox {C}$, and death via thermal stress occurs rapidly by 40 $^{\circ }\hbox {C}$ (A. gambiae survival also tends to increase at higher relative humidity). Several works have relied on a simple survival curve-fit by Martens et al. (1995b) to three data points from 1949 (see Fig. 14), but better data is now available. Bayoh, in a 2001 dissertation (Bayoh 2001), presented data on adult A. gambiae survival in conditions between 5 and 45 $^{\circ }\hbox {C}$, at relative humidities of 40, 60, 80, and 100%, as given in Fig. 8. This data has been used in a number of models, beginning with Mordecai et al. (2013), but these generally assume exponentially distributed survival, which is not actually consistent with the data, as also seen in Fig. 8.

Most recently, Christiansen-Jucht et al. (2014) reared A. gambiae larvae under four different temperatures (23, 27, 31, and 35 $^{\circ }\hbox {C}$), recorded their survival, and then examined the survival curves for adults from these cohorts under the same temperature regimes (no larvae survived at 35 $^{\circ }\hbox {C}$, so this was excluded from adult experiments). Higher temperatures reduced survival of either mosquito stage in general, and adult attrition was exacerbated by a disconnect between larval and adult temperatures. For example, mosquitoes subjected to 31 $^{\circ }\hbox {C}$ as adults were protected somewhat by higher larval temperatures, while adult mosquitoes kept at 23 or 27 $^{\circ }\hbox {C}$ suffered noticeably higher attrition if they emerged from 31 $^{\circ }\hbox {C}$ conditions. Overall, survival was reasonably well-described by a Gompertz distribution, with Gompertz parameters as a function of larval ($T_L$) and adult temperatures ($T_A$), given in Fig. 9.

Fig. 8

Adult A. gambiae survival data adapted from Bayoh (2001). The top left panel gives reported mean female survival as a function of temperature and relative humidity (adapted from Table 6.1 of Bayoh (2001)), while the top right shows this transformed into daily survival probability, assuming exponentially distributed survival times (this was in turn used to inform daily survival in Mordecai et al. (2013)); note that these curves are not greatly dissimilar to Martens’ survival fit (Fig. 14). The bottom left shows survival curves at 80% RH, adapted from Figure 6.4 of Bayoh (2001). The bottom right shows that while a Gompertzian survival function gives an excellent fit to the data (curves shown for 10, 20, and 30 $^{\circ }\hbox {C}$), an exponential distribution fits rather poorly

Fig. 9

Gompertzian survival distributions fit from Christiansen-Jucht et al. (2014), who recorded adult A. gambiae survival at different temperatures following larval maturation at various temperatures

Blood meal dependence Seeking a host and taking a blood meal are much riskier endeavors than resting to digest blood, and therefore Lindsay and Birley (1996) suggested that survival may be relatively constant per gonotrophic cycle (at about 50% per cycle), regardless of cycle length, and Dawes et al. (2009) also observed a spike in mortality after feeding A. stephensi even under experimental conditions. The notion of constant mortality per blood meal was incorporated into a vector lifecycle model by Hoshen and Morse (2004). As a late-night feeding species, A. gambiae may be less likely to suffer attrition when feeding than many other mosquitoes (Klowden and Briegel 1994).

Infection dependence Dawes et al. (2009) observed the survival of A. stephensi-fed Plasmodium-infected blood to decrease with increasing parasite load (as ookinetes), and as reviewed, infection with a higher number of oocytes also seems to incur a survival cost. An older meta-analysis by Ferguson and Read (2002) also concluded that Plasmodium infection decreases mosquito survival, and intriguingly, Pollitt et al. (2013) more recently observed that infection with higher oocyst densities both decreased vector survival and, perhaps via increased competition among parasites for nutrients or a more robust immune response, also decreased the number of infectious sporozoites resulting from infection. In sum, infection is not benign in the adult mosquito, but this has been rarely, if ever, incorporated into malaria transmission models.

5.2.3 Temperature and immature development and survival

Both development time and mortality of aquatic stage anophelines are clearly and nonlinearly related to temperature, although some models (e.g., Craig et al. 1999; Hoshen and Morse 2004; Parham and Michael 2010) have only considered development time as a decreasing function of temperature with mortality temperature-independent. Several works (Craig et al. 1999; Hoshen and Morse 2004; Parham and Michael 2010; Alonso et al. 2011) have used a relation derived from a 1947 work by Jepson et al. (1947), who reported the development rate of A. gambiae as a function of mean temperature at 11 natural breeding sites in Mauritius, with the larval duration time (in days), l(T), reported by Craig et al. (1999) as

$$\begin{aligned} l(T) = \frac{1}{0.00554T - 0.06737}. \end{aligned}$$

(22)

More recently, Bayoh and Lindsay (2003, 2004) determined A. gambiae larval development and survival as a function of temperature under laboratory conditions, at constant temperatures between 10 and 40 $^{\circ }\hbox {C}$. Now, in Bayoh and Lindsay (2003), development time from egg to adult generally decreased with temperature from 18 to 26 $^{\circ }\hbox {C}$, while leveling off at about 10 days from 26 to 32 $^{\circ }\hbox {C}$, as shown in Fig. 10. Below 18 $^{\circ }\hbox {C}$ and above 32 $^{\circ }\hbox {C}$, no eggs survived to adulthood, presumably because of high attrition rates, rather than altered development time. Overall development rate, r(T), was described by the authors using the function

$$\begin{aligned} r(T) = a + b T + c e^T + d e^{-T}. \end{aligned}$$

(23)

In a second work, Bayoh and Lindsay (2004) describe temperature-dependent larval survival time and the fraction surviving to adulthood, as shown in Fig. 10. It is important to note that survival did not follow an exponential distribution. That is, life expectancy was highly correlated with age, with the age-life expectancy curve shifted by ambient temperature, such that for a given temperature, most larvae live a similar lifespan. Therefore, as with adult survival (Sect. 5.2.2), a simple differential equation assuming first-order death kinetics is not a true representation of the biology.

Fig. 10

The left panel gives immature A. gambiae development time as a function of temperature, from Bayoh and Lindsay (2003), and the fit to Eq. (23) ($a = -.05$, $\hbox {b} = 0.005$, $\hbox {c} = -2.139 \times 10^{-16}$; $\hbox {d} = -281357.656$), while the right panel gives overall survival time and the percentage of immature mosquitoes surviving to adulthood (from Bayoh and Lindsay (2004)). Note that the peaks of the survival duration and survival to adulthood curves are offset as a consequence of temperature-dependent development. That is, while there is a penalty to absolute survival duration within the 20–30 $^{\circ }\hbox {C}$ range, immature development time also falls over this temperature range. Thus, up to about 25 $^{\circ }\hbox {C}$, the faster development time outweighs the mortality cost, such that overall survival to adulthood is maximized at higher temperature values than absolute survival time

5.2.4 Larval density and immature development and survival

Larval population density negatively affects anopheline larval survival and other life history traits, primarily development time and adult size, as demonstrated in several laboratory and semi-natural artificial breeding site experiments (Lyimo et al. 1992; Schneider et al. 2000; Gimnig et al. 2002; Jannat and Roitberg 2013; Muriu et al. 2013), although the exact relationships between density and life history vary somewhat with experimental setting. Rainfall and other climate variables strongly determine the size and availability of aquatic habitats, especially the small temporary pools preferred by A. gambiae (Minakawa et al. 1999), and thus, understanding larval density-dependence is necessary for a full accounting of weather and the anopheline lifecycle.

For various mosquitoes, it had generally been found that crowding leads to longer development times, lower survival, and smaller adults, but this was not tested in anophelines until 1992, when Lyimo et al. (1992) raised A. gambiae in the lab in plastic trays at either 24, 27, or 30 $^{\circ }$C, at densities of 0.5, 1.0, or 2.0 larvae/cm$^2$. Yet, their results were rather curious, suggesting increased mortality with increasing density at 24 or 30 $^{\circ }\hbox {C}$, but the opposite at 27 $^{\circ }\hbox {C}$, while higher densities actually seemed to slightly reduce development time. Schneider et al. (2000) performed similar experiments upon populations of A. gambiae and A. arabiensis at 27 $^{\circ }\hbox {C}$, again at densities of 0.5, 1.0, or 2.0 larvae/cm$^2$, and observed the highest density to reduce survival, while the density effect on development time was equivocal.

In sharp contradistinction, three more recent studies (Gimnig et al. 2002; Muriu et al. 2013; Jannat and Roitberg 2013) that used lower larval density ranges showed a very clear negative relationship between density and both development time and adult mosquito weight. Gimnig et al. (2002) raised A. gambiae larvae in outdoor artificial habitats mimicking typical field conditions in west Kenya (essentially dried mud pits), each containing about 1 L of water and 600 cm$^2$ of water surface area, and examined life history parameters across a range of densities, from 0.0333–0.333 larvae/cm$^2$: the development and weight trends were very clear. Furthermore, while survival did decrease with density, this was not significant under statistical analysis. Muriu et al. (2013) performed similar outdoor experiments in coastal Kenya (density range 0.0333–0.5322 larvae/cm$^2$), published density-dependent survival and development curves, and found development rate, survival, and weight to uniformly decrease with density.

Finally, Jannat and Roitberg (2013) most recently attempted to separate the effects of competition for food from crowding per se in A. gambiae, by raising larvae at different densities with food at either a fixed per capita level or at a fixed total level (and at 30 ± 2 $^{\circ }\hbox {C}$, 75–80 %RH). Even with adequate food per larvae, crowding led to higher mortality, smaller adults, and a skewed male:female ratio favoring females. A possible mechanism for the latter is that male larvae are smaller and thus may be more vulnerable to crowding stresses. Under fixed total food resources, time to adulthood was additionally prolonged, and larvae similarly suffered increased mortality and smaller size. Note that crowding, by itself, did not lead to prolonged development time, which may therefore be a consequence solely of nutritional stress.

Given the divergent results of the lower (Gimnig et al. 2002; Muriu et al. 2013; Jannat and Roitberg 2013) and higher density (Lyimo et al. 1992; Schneider et al. 2000) studies, we must then ask, what larval densities are typically encountered in the field? Early field studies performed and reviewed by Service (1971), using 100 mL dippers with a 9.5 cm diameter, gave 0.5–2 fourth-instar larvae per dip. Since about 5% of all larvae were fourth-instar in these studies, this suggests an overall density of 0.035–0.14 larvae/cm$^2$, and is congruent with mark-release-recapture experiments in five pools (Service 1971) that suggested 0.04–0.10 larvae/cm$^2$ (again assuming 5% of larvae are fourth-instar).

Much more recently, Kweka et al. (2012) sampled 51 aquatic habitats in western Kenya over 85 weeks, and found about 6 A. gambaie larvae per 350 mL dipper in hoofprints and swamps. Supposing, say, a circular right angle geometry for hoofprints with a generous depth of 10–25 cm, then this translates into no more than 0.05–0.15 larvae/cm$^2$.

Considering all the aforementioned collectively, it seems likely that, at low to moderately high densities that are within the range typically encountered in the field, the dominant effect of increasing density is nutrient competition, in turn resulting in delayed development, smaller adult sizes, and increased mortality. Developmental effects may plateau around 0.5 larvae/cm$^2$, and at very high densities crowding may directly stress larvae to reduce survival. Survival seems to be affected across density ranges, but the exact relationship varies significantly among studies. Development time as a function of larval density from the studies reviewed above are compared graphically in Fig. 11, along with survival curves adapted from Muriu et al. (2013). In conclusion, those lab/semi-field studies employing lower densities are much more realistic and relevant, and thus we may expect density to deleteriously affect both survival and development rates, although often only the former is considered in models, but see, e.g., Lunde et al. (2013b) for an exception.

Fig. 11

The left panel gives the relationship between larval density (per unit surface area) and mean development time, with data extracted from five studies (Gimnig et al. 2002; Muriu et al. 2013; Jannat and Roitberg 2013; Lyimo et al. 1992; Schneider et al. 2000). Curves are labeled with the approximate experimental mean temperature, and there are three such curves from Lyimo et al. (1992). Note that mean development time for Muriu et al. (2013) was extracted from the survival curves given on the right of the figure, where the top gives cumulative portion of larvae surviving to pupation, and the bottom shows time to pupation for those that survived to pupation

5.2.5 Rainfall: simple models for carrying capacity, oviposition, and survival

Modeling of rainfall’s effect on the vector lifecycle is more variable than that of temperature in published works, and it is more uncertain. By determining the quantity of habitat available, rainfall affects both oviposition and density-dependent larval development and survival. Excessive rainfall also can increase immature mosquito attrition via washout of habitats, and this has been modeled as well. We review several of these phenomena and basic modeling approaches here, while we explore more complex physical modeling of the anopheline microhabitat in the next section.

Simple carrying capacity models Perhaps the most straightforward way to model this notion is to have the carrying capacity of larval habitats, K, be a function of recent rainfall. Yé et al. (2009) modeled growth of the adult mosquito population very simply using a logistic growth term, with K linearly proportional to the prior week’s summed rainfall (under the assumption that several days’ worth of rain contribute to breeding pools); the larval compartment was not modeled. Several other recent works used a ‘quasi-logistic” term to describe the larval death rate (White et al. 2011; Christiansen-Jucht et al. 2015). This was incorporated into a simple mosquito lifecycle model by White et al. (2011) as

$$\begin{aligned} \frac{dE}{dt}= & {} \beta M - \frac{1}{d_E} E - \mu _{E_0} E \left( 1 + \frac{E + L}{K(t)}\right) , \end{aligned}$$

(24)

$$\begin{aligned} \frac{dL}{dt}= & {} \frac{1}{d_E} E - \frac{1}{d_L} L - \mu _{L_0} L \left( 1 + \gamma \frac{E + L}{K(t)}\right) , \end{aligned}$$

(25)

$$\begin{aligned} \frac{dP}{dt}= & {} \frac{1}{d_L} L - \frac{1}{d_P} P - \mu _P P, \end{aligned}$$

(26)

$$\begin{aligned} \frac{dM}{dt}= & {} \frac{1}{2} \frac{P}{d_P} - \mu _M M, \end{aligned}$$

(27)

where E and L are the number of early and late larval instars, respectively, P is pupae count, and M is the adult female mosquito count; $d_i$ is the mean duration of stage i, and the 1 / 2 factor in Eq. (27) accounts for the emerging adult male:female sex ratio. The larval “carrying capacity”, K(t) (in units larvae), was modeled as a convolution of recent rainfall with some weighting function, either a constant, linearly decreasing, or exponentially decreasing function; the latter is given as

$$\begin{aligned} K(t) = \lambda \frac{1}{\tau \left( 1 - \exp \left( \frac{-t}{\tau }\right) \right) } \int _0^t \exp \left( \frac{-(t - \acute{t})}{\tau }\right) \text {rain}(\acute{t}) d\acute{t}, \end{aligned}$$

(28)

where $2\tau$ is the mean of the exponential distribution, $\text {rain}(t)$ is daily rainfall, and $\lambda$ is a free scalar. In other terminology, K(t) is given by passing rainfall through a low-pass filter, with one option (the exponential filter) a leaky integrator.

It is important to point out that, under this formulation, K(t) is not actually a carrying capacity. It is, rather, the inverse of a density-dependent death term, as expanding Eq. (2526), for example, gives

$$\begin{aligned} \frac{dE}{dt} = \beta M - \frac{1}{d_E} E - \mu _{E_0} E - \mu _{E_0} \frac{E^2}{K(t)} - \mu _{E_0} E \frac{L}{K(t)}. \end{aligned}$$

(29)

The actual carrying capacity is a complex function of K(t), $\beta$, and the other model parameters. One can see this should be true by analogy to the simpler Verhulst equation (Vogels et al. 1975) for population growth, which imagines death to increase with the square of population as

$$\begin{aligned} \frac{dy}{dt} = a y - b y^2, \end{aligned}$$

(30)

which can be rearranged to the “ecological” logistic equation as

$$\begin{aligned} \frac{dy}{dt} = r y \left( 1 - \frac{y}{K}\right) , \end{aligned}$$

(31)

with $r = a$ and $K = a/b$ the “carrying capacity”. Thus, carrying capacity here is not independent, but is a function of other more fundamental model parameters. It is important that claims concerning the biological meaning of mathematical terms are internally consistent, otherwise serious errors could be introduced, say, when parameterizing a “carrying capacity” from experimental data or in the interpretation of model output. Several other models (such as Agusto et al. (2015), Okuneye and Gumel (2017)) have used the logistic term to model the birth rate of immature mosquitoes (although without explicit dependence of the carrying capacity on rainfall, but see Depinay et al. (2004) for an exception), a setting where it more properly imposes a well-defined limit to larval population size.

Oviposition Oviposition is affected by both the raw availability of breeding habitat, and the density of larvae within potential positing sites (Sumba et al. 2008). Hoshen and Morse (2004) modeled this very simply, by assuming that each oviposition event yields $\gamma R_d$ eggs, where $R_d$ is the dekadal (i.e., sum of the prior ten days) rainfall in mm, and $\gamma$ was 1 egg/mm. It is reasonable that rainfall over the recent past should be integrated, as it is the sum of standing water that provides habitat, but this simple linear relation is probably not supported.

Sumba et al. (2008) studied experimentally how the presence of larvae in aquatic habitat either encouraged or discouraged A. gambiae oviposition. Interestingly, they found that while pre-existing larvae uniformly discouraged oviposition when distilled water was used, when natural anophelene pool water was used, low densities of larvae actually encouraged oviposition with a shift to deterrence only at high density. Additionally, the larger late instars were more of a deterrent, while the presence of one-day old eggs had no effect either way. As illustrated in Fig. 12, the experimental assay offered gravid mosquitoes a choice between two pools (2 cm deep, 4 cm diameter, 12.57 cm$^2$ SA, 20 mL liquid volume), one containing between 1 and 40 early or late instars (test), the other empty (control). Pool preference was quantified by the oviposition index (OI), defined as

$$\begin{aligned} \text {OI} = \frac{N_t - N_s}{N_t + N_s}, \end{aligned}$$

(32)

where $N_t$ is the number of eggs laid in the test pool and $N_s$ the number in the control pool. The OI ranges from $-1$ to $+1$, with positive values indicating a preference for the test pool, and Sumba et al. (2008) fit oviposition data using the functional form

$$\begin{aligned} \text {OI} = a \exp \left( -b (x-c)^2\right) - d, \end{aligned}$$

(33)

where x is the number of larvae and a, b, and c are free parameters; results for early and late instars are given in Fig. 12. Now, these experiments used extraordinarily high larval densities, as field studies reviewed previously (Service 1971; Kweka et al. 2012) suggest less than 0.02 larvae/mL is more common, and therefore, within plausible field densities larval density likely has no or a slightly positive effect on oviposition.

Fig. 12

Oviposition index in A. gambiae as a function of larval count (in 20 mL of liquid) in the test pool, when larvae are either early (left panel) or late (right panel) instars, as determined by the experiments of Sumba et al. (2008) where water was taken for natural anophelene pools (the presence of larvae uniformly deterred oviposition when distilled water was used). Curves are fits to OI = $a \exp (-b (x - c)^2) - d$, where x is the number of larvae, and for early instars $a = 1.037$, $b = 0.015$, $c = 6.34$, $d = 0.616$; for late instars $a = 1.1136$, $b = 0.00315$, $c = 6.42$, $d = 0.9524$

These results were incorporated into a model by Parham et al. (2012), who recast Sumba et al.’s curve fits for OI to depend upon density, not larval count, with density calculated from total habitat volume, itself determined according to a hydrodynamics model (presented in part below), and then assumed that the fraction of eggs a gravid female lays, $f_t$, is

$$\begin{aligned} f_t = \frac{N_t}{N_s + N_t} = \frac{1}{2}(\text {OI} + 1). \end{aligned}$$

(34)

While this relation has an empirical basis, it is unclear if only half of eggs should be laid at an OI of 0, which simply implies no preference for empty over already occupied pools, or, more generally, that a preference rank for unoccupied pools in a simple forced choice experiment translates linearly to overall oviposition likelihood in a general environment.

Rainfall-dependent mortality A widely-cited 2007 experiment by Paaijmans et al. (2007) placed either first- or fourth-instar A. gambiae larvae in outdoor artificial basin habitats in western Kenya over the course of the rainy season, and observed significantly increased rates of both flushing losses and larvae mortality during rainy nights, with the younger first-instar larvae’s suffering greater, in absolute terms. On this basis, several models have posited different functional forms for increased mortality from rainfall, with Parham and Michael (2010), for example, giving a unimodal relationship between egg survival and rainfall (see Sect. 6.2.2), and Tompkins and Ermert (2013) took larval mortality to increase with precipitation.

5.2.6 Rainfall, hydrodynamics, and the microhabitat

Several more complex, but more realistic, physical descriptions of the microhabitat geometry and heat and water balance have been proposed (Paaijmans et al. 2008a, b; Parham et al. 2012; Asare et al. 2016a, b, c); models for regional-scale hydrology have also been applied to estimating malaria transmission (Bomblies et al. 2008, 2009; Bomblies 2012; Tompkins and Ermert 2013; Asare et al. 2016c), but we restrict our attention here to microscale dynamics. Detailed microscale models have multiple advantages over the simpler approaches discussed above. First, rainfall can directly inform habitat volume and surface area, yielding a time-dependent immature carrying capacity (broadly defined) from basic physical principles and geometric parameters. Second, such models relate water temperature to ambient air temperature in a physically realistic and non-constant manner. Finally, one may generate predictions on how local variations in habitat geometry, such as shading, interact with more global parameters, such as ambient temperature. In this section, we review the basic construction of a comprehensive microhabitat model. Figure 13 summarizes the key parameters describing habitat geometry, and the major mechanisms for heat and water volume loss/gain.

Before continuing, we also note that for any hydrodynamic model describing habitat volume and/or surface area, such metrics must be translated into some kind of carrying capacity or density-dependent death term, etc. in the mosquito population dynamics model, and several authors have assumed a biomass carrying capacity for anopheline ponds to be about 300 mg m$^{-2}$, with fourth stage instars weighing 0.45 mg (Depinay et al. 2004; Bomblies et al. 2008; Tompkins and Ermert 2013). A separate method is that of Lunde et al. (2013b), who calculated an immature anopheline carrying capacity at a relatively high spatial scale as a composite function of soil moisture and potential river length, with the latter determined from the HydroSHEDS database, which in turn gives water accumulation potential based upon the Earth’s topology.

Fig. 13

General schematic for a single Anopheles breeding site microhabitat, and the coupled mechanisms determining the overall heat and water balance

Basic geometry To model an Anopheles microhabitat, one must prescribe some volume (V)-area (A)-depth (h) relationship, with one popular option a simple set of equations developed by Hayashi and colleagues (Hayashi and Van der Kamp 2000; Brooks and Hayashi 2002), that describe small topographic depressions in terms of three empirical parameters: maximum surface area $A_{max}$, maximum depth, $h_{max}$, and a dimensionless shape parameter p, such that $p < 1$ and $p > 1$ indicate concave and convex geometries, respectively. If depth, h, is known, A and V are determined as

$$\begin{aligned} A= & {} A_{max} \left( \frac{h}{h_{max}}\right) ^{\frac{2}{p}} \end{aligned}$$

(35)

$$\begin{aligned} V= & {} \frac{A_{max} h_{max}}{1 + \frac{2}{p}}\left( \frac{h}{h_{max}}\right) ^{1 + \frac{2}{p}}. \end{aligned}$$

(36)

If volume, V, is prescribed instead, it is straightforward to rearrange the equations to solve for h and A. Alternatively, Parham et al. (2012) employed a simple right-angle cone as a microhabitat geometry.

Heat-balance The heat and water volume balance within a microhabitat are linked, and, generally speaking, we have the change in heat, dQ / dt (in W), as

$$\begin{aligned} \frac{dQ}{dt} = A (R_n - \lambda E - H - G) + P_Q - I_Q \end{aligned}$$

(37)

where $R_n$ is net radiation per unit area (W m$^{-2}$), $\lambda E$ is latent heat flux (W m$^{-2}$), H is sensible heat flux (W m$^{-2}$), G is heat flux through the surrounding soil (W m$^{-2}$) (Allen et al. 1998; Asare et al. 2016b), while heat is also contained in the water gained via precipitation and runoff, $P_Q$, and lost via infiltration, $I_Q$. We briefly examine each component of Eq. (37). First, net radiation, $R_n$, decomposes as the sum of incoming shortwave radiation, $R_s$, and incoming and outgoing longwave radiation, $L_{in}$ and $L_{out}$, respectively,

$$\begin{aligned} R_n = R_s - L_{out} + L_{in}, \end{aligned}$$

(38)

Shortwave radiation is the fraction of the global horizontal irradiance (GHI) that is not reflected or blocked by shade,

$$\begin{aligned} R_s = \text {GHI} (1 - a) (1 - SF), \end{aligned}$$

(39)

where a is the albedo of water and SF is a shade factor. The longwave radiation balance can be determined from the relatively simple relations,

$$\begin{aligned} L_{out} = \epsilon _w \sigma T_{w(K)}^4, \end{aligned}$$

(40)

$$\begin{aligned} L_{in} = \epsilon _a \sigma T_{a(K)}^4 (1 - SF), \end{aligned}$$

(41)

where $\epsilon _w \approx 0.98$ the emissivity of water (Paaijmans et al. 2008b), $\sigma$ = $5.67 \times 10^{8}$ W m$^{-2}$ K$^{-4}$ the Stefan-Boltzmann constant, $T_{w(K)}$ and $T_{a(K)}$ are water and air temperatures in Kelvin, and finally, $\epsilon _a$ is either the clear-sky emissivity, or the cloud-corrected atmospheric emissivity. A variety of algorithms exist to estimate $\epsilon _a$, as reviewed by Flerchinger et al. (2009), with one simple option for clear-sky emissivity due to Ångström (Flerchinger et al. 2009) given as

$$\begin{aligned} \epsilon _{a} = \left( .83 - .18 \times 10^{-.067 e_a}\right) , \end{aligned}$$

(42)

where $e_a$ is saturation pressure (kPa).

Latent heat flux, $\lambda E$, determined from the mass of water evaporated per unit time, E (kg day$^{-1}$), and the latent heat of evaporation, $\lambda$, equal to 2.45 MJ kg$^{-1}$ at 20 $^{\circ }\hbox {C}$ (Allen et al. 1998), is a relatively complex phenomenon, as it involves both storage of heat into the latent form of water vapor at the water surface, and the removal of this vapor from the surface. In general, with lower vapor pressure at the water surface and faster wind speed, both processes (water vapor formation and removal) are accelerated. These qualitative notions can be simply formalized to model evaporation (i.e. water mass loss) as a bulk transfer process (Sene et al. 1991; Paaijmans et al. 2008a; Asare et al. 2016b)

$$\begin{aligned} E = C u (e_{sw} - e_a), \end{aligned}$$

(43)

where C is the mass transfer coefficient, u (m s$^{-1}$) is wind speed at reference height, $e_{sw}$ is the vapor pressure at saturation for the water surface temperature, and $e_a$ is the atmospheric vapor pressure at reference height. An alternative method is to use the Penman-Monteith equation, or a variation thereof, derived from the simultaneous solution of equations for energy balance and mass transfer as (Finch and Hall 2001)

$$\begin{aligned} E = \frac{1}{\lambda } \left( \frac{\varDelta \varLambda + \frac{\rho c_p (e_s - e_a)}{r_a}}{\varDelta + \gamma \left( 1 + \frac{r_s}{r_a}\right) }\right) \end{aligned}$$

(44)

where $\varLambda$ is the available energy (typically taken as $R_n - G$ (Allen et al. 1998)), $\varDelta$ is slope of the vapor pressure curve (kPa $^{\circ }\hbox {C}$ $^{-1}$), $e_s$ is the saturation vapor pressure, $e_a$ is the actual vapor pressure, $\rho _a$ is the density of air, $c_p$ is the specific heat of air, $\gamma$ is the psychrometric constant, defined as

$$\begin{aligned} \gamma = \frac{c_p P}{\epsilon \lambda }, \end{aligned}$$

(45)

where P is the atmospheric pressure and $\epsilon$ = 0.622 is the ratio of the molecular weight of water vapor to dry air. Furthermore, in this formulation, water vapor mass transfer is governed by two resistances in series, first a “bulk” surface resistance, $r_s$, which is a sum measure of the resistance to vapor flow from soil and vegetation and may be set to zero over open water, and a second aerodynamic resistance, $r_a$, which is inversely proportional to wind speed (see Finch and Hall (2001) or Allen et al. (1998)). Parham et al. (2012) employed the FAO, or modified, Penman-Monteith equation for transevaporative flux from a vegetated surface to describe water loss from anopheline habitat; further details may be found in Parham et al. (2012) and Allen et al. (1998).

Similar to Eq. (43), sensible heat flux, H, may be given as (Asare et al. 2016b)

$$\begin{aligned} H = \rho _a c_p C u (T_w - T_a), \end{aligned}$$

(46)

and G is typically taken as some small fraction, f, of $R_n$, perhaps 0.15 (Asare et al. 2016b; Paaijmans et al. 2008a), i.e.

$$\begin{aligned} G= f R_n. \end{aligned}$$

(47)

Finally, the heat contained in precipitation and infiltration water ($P_Q$ and $I_Q$ respectively) is simply determined from the density and specific heat of water, and using the volume-balance principles presented next.

Volume-balance The change in water volume, dV / dt, may be approximated as a function of an imposed precipitation time-series, P(t) (m day$^{-1}$ or m s$^{-1}$ if appropriate), and loss to evaporation (E) and infiltration (I), as

$$\begin{aligned} \frac{dV}{dt} = P(t) (A + R_{frac} (A_{catch} - A)) - A (E + I), \end{aligned}$$

(48)

where $A_{catch}$ is the catchment area for precipitation runoff, $R_{frac}$ is that fraction of runoff water within the catchment area that makes it to the habitat, and Eq. (48) is also subject to the constraint that V not exceed $V_{max}$. Evaporation, E, is determined as above. Infiltration of water in sandy pools in the Sahel region proceeds in a roughly biphasic manner, where water is initially lost rapidly to the porous sandy soil, but low permeability clay that collects at the bottom creates a “clogged” zone, in which infiltration slows dramatically (Desconnets et al. 1997; Porphyre et al. 2005). This phenomenon was simply modeled by Asare et al. (2016a), who gave I as

$$\begin{aligned} I = I_{max} \left( \frac{A}{A_{max}}\right) , \end{aligned}$$

(49)

where $I_{max}$ (m day$^{-1}$) is the maximum infiltration rate. A more comprehensive model could also incorporate overflow and washout of immature anophelines during heavy rains.

Water versus air temperature Water and air temperatures in open anopheline habitats are likely to vary by 3–6 $^{\circ }\hbox {C}$ or more (Parham et al. 2012), and such a disparity could strongly affect optimal ambient temperatures for malaria transmission. Paaijmans et al. (2008b) observed mean water temperatures to be several degrees Celsius above the surrounding air in artificial habitats, with the difference greatest during the day, when solar gain into pools is high. The difference between maximum air and water temperatures could exceed 10 $^{\circ }\hbox {C}$, and, unsurprisingly, there was greater thermal stability in larger pools (see also Paaijmans et al. (2008a)). Such an air-water temperature disparity is also observed in simulations of the theoretical microhabitat framework presented above. Thus, while many works assume a constant air-water temperature difference (often zero), this is probably overly simplistic.

6 A partial genealogy of recent weather-driven malaria models

It is well beyond our scope to review all malaria mathematical models that incorporate weather, and we must omit any discussion of many excellent works (a partial list of works not considered further here includes (Depinay et al. 2004; Lou and Zhao 2010; Cailly et al. 2012; Dembele et al. 2009; Bomblies et al. 2009; Eckhoff 2011; White et al. 2011; Lunde et al. 2013a, b; Tompkins and Ermert 2013; Nikolov et al. 2016)). We restrict our attention to process-based, mechanistic models, and among these focus on several influential lines of work from the past two decades, beginning with widely cited works from the 1990s by Martens and colleagues which suggested a significantly increased malaria range with global warming, and moving through a 2013 paper by Mordecai et al. (2013), who concluded prior works had significantly overstated the optimum temperature range for transmission. Beyond that work, we highlight several recent efforts with slightly different focuses, namely immunity (Agusto et al. 2015; Yamana et al. 2013, 2017), host age-structure (Okuneye and Gumel 2017), and age-dependent vector survival (Christiansen-Jucht et al. 2015). Recently, the importance of daily and seasonal temperature fluctuations in determining malaria potential has recognized, and we discuss several recent works which make this their focus (Paaijmans et al. 2009, 2010, 2013b; Blanford et al. 2013; Beck-Johnson et al. 2017). All models are informed by the empirical vector/parasite-weather relations covered in the prior section, and many directly employ the Ross–Macdonald model framework. We close this section with a discussion of host mobility in multi-patch geometries, which heretofore has not been incorporated in weather-driven malaria models, but is of potentially great value in studying the possible combined roles of climate change and human mobility in the spread of malaria into highland Kenya.

6.1 Epidemic potential models (1995–1999)

Several works by Martens and others (Martens et al. 1995a, b, 1997, 1999), published in the late 1990s, used the notion of “epidemic potential” (EP), defined to be the inverse of the critical vector density, $m_c$ (units mosquitoes/man), in turn derived from the Macdonald model by setting ${\mathcal R}_0 = 1$ (see Eq. (15)) and solving for m, giving

$$\begin{aligned} m_c = k_1 \frac{-\ln (p)}{a^2 p^n}, \text {EP} = \frac{1}{m_c}, \end{aligned}$$

(50)

where $k_1$ is a constant, and equal to r / (bc) for Macdonald’s model. The parameters n (sporogonic duration), a (daily biting rate), and p (daily mosquito survival) were determined as a function of temperature, with n determined via Moshkovsky’s formula ($D = 111$ and $T_{min} = 16-19$), a assumed to be proportional to the length of the gonotrophic cycle and similarly determined from Moshkovsky’s formula (Eq. (18), with $D = 36.5$ and $T_{min} = 9.9$), and p determined via a trinomial fit to three data points from a 1949 work, such that

$$\begin{aligned} p(T) = \exp \left( \frac{-1}{-4.4 + 1.31 T - .03 T^2}\right) . \end{aligned}$$

(51)

Fig. 14

Temperature dependence of vector and parasite parameters as per Martens et al. (1995b), Martens et al. (1997). The top panels show the duration of the sporogonic and gonotrophic cycles as functions of mean temperature, according to the formula of Moshkovsky using $D = 111$ and $T_{min} = 16$ and $D = 36.5$ and $T_{min} = 9.9$, respectively. The bottom shows adult mosquito survival per the relation derived by Martens et al. (1995b), and the normalized epidemic potential (see Eq. (50)), assuming the inverse of the biting rate (1 / a) is equal to the gonotrophic duration

These parameters, and the resulting EP, are graphed as functions of temperature in Fig. 14. Using geographical data from GCMs, these authors then concluded that global warming could cause overall epidemic potential to increase by 12–27% in 2050 (Martens et al. 1997). The idea of EP is primarily applicable to areas not currently endemic for malaria: it is meant to elucidate what regions may become vulnerable in the future. However, Rogers and Randolph (2000) pointed out that, in areas where ${\mathcal R}_0 < 1$ (and, thus, not vulnerable to epidemics), an increase in EP does not mean that ${\mathcal R}_0$ increases above 1. Thus, the EP notion can overestimate the effect of temperature changes. Even considering that critical mosquito densities may decrease with warming, this is a very incomplete picture, as it merely gives a threshold mosquito density necessary to spark an epidemic, while failing to predict how mosquito densities will change under climate change. This is key: we cannot assume that altering major anopheline lifecycle parameters will affect disease transmission but not the Anopheles population itself.

The relations reported by Martens et al. (1995b) helped inform a widely cited effort by Craig et al. (1999) (see also Snow et al. (1999)) that used “fuzzy logic” to derive maps of climatic suitability for malaria transmission in Africa. These authors determined that, when rainfall is not limiting, yearly mean temperatures above 22 $^{\circ }\hbox {C}$ lead to perennial infection, 18 $^{\circ }\hbox {C}$ is too cold for stable transmission but does allow epidemics in warmer years, while 15 $^{\circ }\hbox {C}$ is prohibitory.

A later 2011 work by Gething et al. (2011) also used similar temperature-dependent terms and the notion of vector capacity to map global temperature constraints on P.falciparum and P.vivax transmission.

6.2 More complex models through Mordecai et al. (2013)

6.2.1 Hoshen and Morse model (2004), reformulation and extensions (2011)

In 2004, Hoshen and Morse (2004) developed a more comprehensive model of the mosquito lifecycle, explicitly including temperature dependent progression through the immature (egg, larvae, pupae) mosquito stages (considered as one), coupled with the temperature-dependent gonotrophic cycle of adult female mosquitoes, by which a blood meal is taken to yield oviposition of new eggs (the number of eggs laid is dependent upon the sum of the prior ten day’s rainfall), along with a sporogonic cycle and a basic susceptible-exposed-infected-susceptible (SEIS) model for human infection. Note that the sporogonic cycle in infected mosquitoes advances independently of the gonotrophic cycle, and the overall model architecture is schematized in Fig. 15. Immature mosquitoes progress through physiologic time (as opposed to chronological time) at a rate, m(T) (1/day), determined as the inverse of the sum of the duration of all larval stages (note that the immature mosquito class is still considered as one by Hoshen and Morse (2004)), with these durations determined from Jepson et al. (1947) (Eq. (22) of Sect. 5.2.3), which gives a hyperbolic relationship between total time in the immature stage and temperature. Death occurs daily, independent of temperature, with an assumed daily survival of 90%. The model follows previous work in using Moshkovsky’s hyperbolic formula for sporogonic, $S_C$, and gonotrophic, $G_C$, durations (see Eq. (18) and Sects. 5.1.1 and 5.2.1).

Fig. 15

A schematic for the essential elements of the Hoshen and Morse model. Note that this schematic represents our particular interpretation of the original difference equation-based model

At the end of each gonotrophic cycle, the (adult female) mosquito takes a blood meal, surviving the risky adventure with probability $\alpha \approx 0.5$, and thus daily survival probability is $\alpha ^{1/G_C}$, giving an overall per-capita death rate of $-\ln (\alpha ^{1/G_C})$ (denoted $\mu$ in the model equations below). Note that adult survival is therefore indirectly coupled to temperature, and no other mechanism for death is considered. After each gonotrophic cycle, mosquitoes lay $\gamma$ $R_d$ eggs, where $R_d$ is the sum of the prior ten days of rainfall in mm, and $\gamma$ was 1 egg/mm, thus coupling oviposition to habitat created by rain. The fraction of blood meals taken on humans is B, and mosquitoes contract malaria from infected humans with probability $\chi$. If a mosquito is infected, then the sporogonic cycle initiates, taking 111 degree days.

Finally, the human component is a simple susceptible-exposed-infectious-susceptible (SEIS) model with a 14 day delay from exposure to infectivity, with no superinfection or immunity included. A small influx of infected individuals is included, presumably from migration. A small, constant influx of mosquitoes is also included in the model. The model was formulated by Hoshen and Morse as a fairly extensive difference equation, which we do not repeat here, but the general model framework can be translated into continuous time, a more mathematically popular setting. Now, there are three major continuous-time forms common among malaria models: ordinary-differential equations (ODEs), delay-differential equations (DDEs), and age-structured advection partial-differential equations (PDEs), where the latter describes the movement of populations through infinite-dimensional time and age (age is analogous to space in a traditional advection equation). The simpler ODE version, which is not a “direct” translation (it is not generally possible to directly translate from a continuous age-structured model to an ODE model of finite dimension) can be written as

$$\begin{aligned} \frac{dL}{dt}= & {} \gamma R_d \frac{1}{G_C(T)} (S_M + E_M + I_M) - m(T) L - \sigma L, \end{aligned}$$

(52)

$$\begin{aligned} \frac{dS_M}{dt}= & {} m(T) L + \lambda _M - \chi B \frac{1}{G_C(T)} S_M \frac{I_H}{N} - \mu S_M, \end{aligned}$$

(53)

$$\begin{aligned} \frac{dE_M}{dt}= & {} \chi B \frac{1}{G_C(T)} S_M \frac{I_H}{N} - \frac{1}{S_C(T)} E_M - \mu E_M\end{aligned}$$

(54)

$$\begin{aligned} \frac{dI_M}{dt}= & {} \frac{1}{S_C(T)} E_M - \mu I_M,\end{aligned}$$

(55)

$$\begin{aligned} \frac{dS_H}{dt}= & {} -\,\beta B I_M \frac{1}{G_C(T)} \frac{S_H}{N} + \rho I_H,\end{aligned}$$

(56)

$$\begin{aligned} \frac{dE_H}{dt}= & {} \beta B I_M \frac{1}{G_C(T)} \frac{S_H}{N} - \frac{1}{14} E_H,\end{aligned}$$

(57)

$$\begin{aligned} \frac{dI_H}{dt}= & {} \frac{1}{14} E_H - \rho I_H + \lambda _H, \end{aligned}$$

(58)

where

$$\begin{aligned} \mu = -\ln \left( \alpha ^{1/G_C(T)}\right) , \end{aligned}$$

(59)

and where L is the lumped larval stage; $S_M$, $E_M$, and $I_M$ are the susceptible, exposed, and infectious adult mosquito populations; $S_H$, $E_H$, and $I_H$ are similarly susceptible, exposed, and infectious host populations,with $N = I_H + E_H + I_H$. Latency from exposure to infection is a fixed 14 days in humans, $\beta$ is the probably an infectious bite infects, $\rho$ is the rate at which infected recover (with no immunity), $\lambda _M$ and $\lambda _H$ are respectively the influxes of adult mosquitoes and infected humans from migration, and it assumed that there is no mortality in the human population. Other parameters are as above. Now, to make this translation from Hoshen and Morse’s discretely age-structured model, we have assumed exponentially waiting times in every stage, but as intimated throughout Sect. 5, and addressed more directly with a model by Christiansen-Jucht et al. (2015), this is not generally a valid assumption.

Staying truer to the original work, it is also possible to write the model as a set of either age-structured advection equations, or, equivalently, a delay-differential equation, with that caveat that it necessary to cast the model in physiologic, rather than chronological, time. The translation from ODE to either of these settings is relatively straightforward, but requires care with the boundary conditions.

At this point, with the Hoshen and Morse work, we have a fairly realistic model for the mosquito lifecycle, explicitly including dependence upon blood meals for reproduction, and monotonically increasing relationships between temperature and larval development, sporogony, gonotrophy (and, hence, biting rate), while daily mosquito survival decreases monotonically (as a consequence of mortality occuring during biting, which increases with temperature).

The model of Hoshen and Morse, also referred to as the Liverpool Malaria Model (LMM), was revised modestly by Ermert and colleagues (Ermert et al. 2011a, b), who also performed a more thorough literature review for parameter values, and also applied the model to West African field data (Ermert et al. 2011b). While Hoshen and Morse present their model as being valid only for predicting zones of epidemic malaria (e.g., mesoendemic or hypoendemic area), given that it does not account for immunity, Ermert et al. (2011a) argue that their updated version can apply to endemic zones as well, but the original authors’ claim seems biologically well-founded. In any case, it is extremely likely that without explicitly incorporating the effects of immunity in a robust manner, predictions concerning climate change or other interventions on malaria in endemic zones will be suspect, given the fundamental importance of the immune dynamic to the history and epidemiology of the disease (see Sects. 7.1 and 3.3).

6.2.2 Parham and Michael model (2010)

The next major effort we discuss is that of Parham and Michael (2010), a delay differential equation that couples an SEI model for mosquito dynamics with an SIR model for the human population, and directly adopts the temperature-dependent relations for adult mosquito survival, sporogonic duration, and gonotrophic duration from Martens et al. (1995b) (Sect. 6.1), the larval maturity relation from Jepson et al. (1947) (Eq. (22)), and finally adds an assumed nonlinear relation between rainfall and daily egg survival probability, $p_E(R)$,

$$\begin{aligned} p_E(R) = \frac{4 p_{ME}}{R_{LE}^2} R (R_{LE} - R), \end{aligned}$$

(60)

where $R_{LE}$ is the threshold beyond which no eggs survive due to washout, and $p_{ME}$ is the maximum daily survival fraction. The model governing equations for susceptible, $S_M$, exposed, $E_M$, and infectious, $I_M$ mosquito populations, and similarly named human populations (including recovered, $R_H$, and total human population, N), with time-dependence suppressed except for delay-terms, is given by

$$\begin{aligned} \frac{dS_M}{dt}= & {} \lambda (R,T) - a(T) b_1 S_M \frac{I_H}{N} - \mu (T) S_M, \end{aligned}$$

(61)

$$\begin{aligned} \frac{dE_M}{dt}= & {} a(T) b_1 S_M \frac{I_H}{N} - \mu (T) E_M - a(T) b_1 S_M(t-\tau _M(T)) l_M(T) \frac{I_H(t-\tau _M(T))}{N},\end{aligned}$$

(62)

$$\begin{aligned} \frac{dI_M}{dt}= & {} a(T) b_1 S_M(t-\tau _M(T)) l_M(T) \frac{I_H(t-\tau _M(T))}{N} - \mu (T) I_M,\end{aligned}$$

(63)

$$\begin{aligned} \frac{dS_H}{dt}= & {} -\,a(T) b_2 I_M \frac{S_H}{N},\end{aligned}$$

(64)

$$\begin{aligned} \frac{dI_H}{dt}= & {} a(T) b_2 I_M \frac{S_H}{N} - \gamma I_H,\end{aligned}$$

(65)

$$\begin{aligned} \frac{dR_H}{dt}= & {} \gamma I_H. \end{aligned}$$

(66)

The adult mosquito influx term, $\lambda (R,T)$, is given as

$$\begin{aligned} \lambda (R,T) = \frac{B p_E(R) p_L(T) p_P(R)}{\tau _E + \tau _L(T) + \tau _P}, \end{aligned}$$

(67)

with B here the number of eggs per oviposition per adult, $p_i$ the daily survival of immature stage i (i = E, L, and P for eggs, larvae, and pupae, respectively), and $\tau _i$ the development time for stage i (note that while $\lambda (R,T)$ is independent of the total adult mosquito population in this model, these parameters should more generally be coupled for maximum fidelity to the underlying biology); also note that the governing equation for $E_M$ may technically be omitted, as it is uncoupled from the rest of the model. Rainfall affects egg survival, $p_E(R)$, per Eq. (60), and larval development hyperbolically per Eq. (22) (Jepson et al. 1947). As above, other temperature-dependent parameters assume the relations of Martens et al. (1995b) given in Sect. 6.1.

The overall model framework is thus similar to the ODE version of Hoshen and Morse’s model, as related in Eqs. (525354555657)–(58), but it incorporates rainfall differently, does not directly couple the gonotrophic cycle to oviposition or death, and does not explicitly consider the immature mosquito life stage. Note that, while biting is temperature-dependent, oviposition is not coupled to biting directly.

Parham and Michael (2010) also take recovery refractory to further infection in humans to be an absorbing state, which is unrealistic for malaria in general, but reasonable as an approximation for a single epidemic in a previously unexposed population. Thus, the model is an extension of the ideas contained in the simpler works of Martens et al. and Craig et al. (1999) to design a continuous, delay differential equation SEIR framework that is more amenable to explicit (rigorous) analysis than the original difference equation formulation of Hoshen and Morse (2004), but one that has some minor mathematical infidelities to the underlying biology. Parham and Michael (2010) derived several complex expressions for ${\mathcal R}_0$ under different simplifying assumptions, and arrived at their key prediction: malaria transmission is maximized, both in endemic and epidemic areas, in the 32–33 $^{\circ }\hbox {C}$ temperature range.

These authors and several colleagues also published a 2012 work (Parham et al. 2012) that considered the vector lifecycle in detail, developing a hydrodynamics model relating rainfall to habitat volume. This allows a calculation of immature mosquito density, and hence several detailed expressions for density-dependent oviposition and density-dependent larval mortality.

6.2.3 Alonso et al. model (2011)

Alonso et al. (2011) developed another ordinary differential equations based model for a highland tea plantation in Kenya. Larval maturity rate and daily mosquito survival follow the monotonic Jepson (Eq. (22)) and Martens et al. relations (Sect. 6.1), respectively, as in the prior models, with sporogonic duration also hyperbolically related to temperature. The model was novel compared to earlier works for including symptomatic and asymptomatic carriers of infection, superinfection of asymptomatic carriers to the symptomatic state, and clinical treatment of symptomatic sufferers. These authors also introduce a new aquatic stage death term, given as the sum of a constant background, temperature-dependent, and rainfall dependent term:

$$\begin{aligned} \delta _L = \delta _0 + \delta _L(T) + \delta _L(R). \end{aligned}$$

(68)

For temperature dependence, $\delta _L(T)$, Alonso et al. (2011) fit a rather complex piecewise defined function to much more modern data than previously used, namely laboratory A. gambiae survival under different temperatures, by Bayoh and Lindsay (2004). However, a fourth-order polynomial also provides a satisfactory fit, and this data is reviewed in Sect. 5.2.3. Rainfall dependence is given by the assumed function

$$\begin{aligned} \delta _L(R) = \delta _R \varTheta (R - \langle R \rangle _{12}), \end{aligned}$$

(69)

where $\delta _R$ is a constant, $\varTheta (x)$ is x if $x>0$ and 0 else, and $\langle R \rangle _{12}$ is a 12-month rainfall moving average. This function is intended to capture washout mortality from heavy rains. As in other works, the biting rate is assumed to be equal to the inverse of the gonotrophic duration, which, in a departure from the Moshkovsky formula (Eq. (18)), is determined from data collected by Afrane et al. (2005) as

$$\begin{aligned} t_a = \frac{1}{0.091678 T - 1.7982}, \end{aligned}$$

(70)

although this preserves the same basic hyperbolic relationship between temperature and gonotrophy. Finally, it is assumed that each oviposition event yields 66 eggs on average.

The application of this model is rather novel, compared to much of the literature, in that it has as its relatively narrow focus a malarious tea plantation under fairly constant conditions over multiple decades, which experienced rising temperatures since the 1970s, and for which detailed time-series data of clinical malaria cases was available dating from that time. The model was trained using data up to 1985, and was then used to generate counterfactual time-series for more recent malaria burden (with confidence intervals), one where the observed warming did occur, and another where temperature in the 1990s remained “similar” to the 1970s. The model clearly predicted that actually observed cases (which were in the range projected with warming), would have been very unlikely without warming, supporting a clear role for temperature in increasing malaria burden in at least one specific area.

6.3 Some recent works (2013–2017)

6.3.1 Beck-Johnson et al. (2013) model

Beck-Johnson et al. (2013) also developed a delay-differential model for immature mosquito development through egg, larva, and pupa, with development time inversely proportional to temperature (according to a power law, based on Lardeux et al. (2008)), but with survival at each immature (and adult) stage, determined according to a Gaussian function of temperature, such that survival is maximized at intermediate temperatures. The mathematical difficulties of non-constant delays in development of each stage were resolved by re-scaling the model to physiologic, instead of chronological, time, and density-dependent larval mortality was also included.

Anticipating Mordecai et al. (2013) somewhat, the model suggested that adult mosquito prevalence is maximized around 20–30 $^{\circ }\hbox {C}$, while the potentially infectious mosquito population (i.e. that fraction of the adult population that survives long enough for at least one temperature-dependent sporogonic cycle to elapse), is maximized at a slightly higher range (24–30 $^{\circ }\hbox {C}$), but both populations abruptly decline beyond 32 $^{\circ }\hbox {C}$.

6.3.2 The Mordecai et al. response (2013)

In 2013, Mordecai et al. (2013) aggregated an updated collection of data to derive functions relating vector and parasite parameters to temperature that were uniformly unimodal (i.e. with a peak at some intermediate temperate), rather than monotonic, as in most of the previous works discussed. Based on this, and using a formula for ${\mathcal R}_0$ derived partly from Parham and Michael (2010),

$$\begin{aligned} \mathcal {R}_0 = \left( \frac{a(T)^2 b c(T) \exp (-\mu (T)/PDR(T)) EFD(T) p_{EA}(T) MDR(T)}{N r \mu (T)^3}\right) ^{\frac{1}{2}}, \end{aligned}$$

(71)

where N is human density, $\mu (T)$ is the adult mosquito death rate, PDR(T) is the sporogonic rate, MDR(T) is the larval development rate, $p_{EA}(T)$ is probability of survival to adulthood, a(T) is the biting rate, and bc(T) is vector competence, these authors concluded that previous works had dramatically overestimated the optimum temperature range for malaria transmission, concluding that transmission is maximized at 25 $^{\circ }\hbox {C}$.

It is important to examine in some detail the thermal response functions used by Mordecai et al. (2013) and their sources. All are graphically illustrated in Fig. 16, with data sources summarized in the caption. Functions are given as either a quadratic polynomial, or according to a three-parameter (c, $T_m$, $T_0$) unimodal function developed by Briere et al. (1999) to describe arthropod development rates:

$$\begin{aligned} r(T) = c T (T - T_0) (T_m - T)^{\frac{1}{2}}. \end{aligned}$$

(72)

Fig. 16

Unimodal thermal-response curves for vector and parasite parameters used by Mordecai et al. (2013); parameters also used by Martens et al. (1997) are also shown in gray. Egg-to-adult survivorship and development rates were fit to Bayoh and Lindsay (2003), adult mortality was estimated from Bayoh (2001), assuming exponentially distributed survival times, and bite rate was estimated from gonotrophic cycle duration data by Lardeux et al. (2008). Vector competence, $b \times c$, (which should be < 1) was estimated from a 1940 work. The rate of the sporogonic cycle was also estimated from several older works, and a newer work by Eling et al. (2001). Fecundity (eggs/female/day) was determined from a work on Aedes albopictus (Delatte et al. 2009), a dengue vector, and is obviously not coupled to the biting rate

One problem with collating these thermal-response curves from different data sources is that the parameters considered are not necessarily independent. For example, larval development is supposed to cease by 34 $^{\circ }\hbox {C}$, yet this may represent a conflation of mortality with development: obviously no larvae complete development if temperatures prohibit survival, but the development rate per se does not necessarily go to zero. As survival and larval development compound multiplicatively in Mordecai et al.’s expression for ${\mathcal R}_0$, this could bias the optimal temperature range downward. A similar argument applies to sporogonic cycle duration, and several other parameters, such as biting rate and eggs laid per day, are also not independent. This basic problem is not unique to the 2013 Mordecai et al. work, and applies deeply to the foundational Ross–Macdonald models.

Lunde et al. (2013a), using an ODE model for malaria transmission, compared six different temperature-dependent adult mosquito death rates, including that of Martens et al. (1997), Ermert et al. (2011a), Parham et al. (2012), Mordecai et al. (2013), and a separate work by Lunde et al. (2013b), all of which, except those due to Martens et al., were determined from data by Bayoh (2001), and all, except an early polynomial also due to Martens et al., similarly suggested transmission to be optimal in the 24–27 $^{\circ }\hbox {C}$ range.

6.3.3 Ryan et al. (2015): Malaria mapping under Mordecai’s curves

In an interesting follow-up, Ryan et al. (2015b) developed a series of malaria potential maps across Africa under current conditions and under projected warming through 2080 (under a mid-range emissions scenario, SRES A1), using the thermal-response functions and ${\mathcal R}_0$ expression of Mordecai et al. (2013) (Equation 71). This work was novel in that it considered both year-round and seasonal malaria potential, and considered the populations, not just land area, at risk. Furthermore, although this work focused on temperature only, while many prior similar mapping efforts have filtered results geographically by applying a minimum rainfall threshold for malaria transmission (e.g. Craig et al. (1999)), Ryan et al. (2015b) applied an aridity mask on the basis of the normalized Normalized Difference Vegetation Index (NDVI), an index determined from satellite measurements of the radiation reflected by a surface, and calculated from the difference in reflectance in the visible light and near-infrared spectrum bands. Values above 0.2 quantify vegetation greenness (Baeza et al. 2011), and at mid- and low latitudes, NDVI correlates with wetness (Suzuki et al. 2006); NDVI has an advantage over rainfall in that it may capture areas of low rainfall still suitable for anophelines due to irrigation, rivers, or other permanent water sources (Baeza et al. 2011; Ryan et al. 2015b).

Ryan et al. (2015b) predicted that, as the globe warms, the areas most suitable for intense, year-round transmission will shift southeasterly from western coastal Africa, with the new peak in transmission potential importantly centered around heavily populated areas such as western Uganda, northern Tanzania, the Lake Victoria region near the Kenyan highlands, and highland Madagascar. They also concluded that, while the total area at any risk for malaria may increase slightly, the area at the highest risk will drop. Of note, this work did not consider daily temperature variations (see Gething et al. (2011) for an example of a mapping effort which does this), and did not consider more detailed hydrodynamics.

6.3.4 Synthesis models with partial immunity and age-structure

Recently, Agusto et al. (2015) proposed a fairly comprehensive ODE model that includes the larval vector stage, with density-limited growth (via a logistic growth term), adapts the temperature-dependent lifecycle parameters from Mordecai et al. (2013), and extended a basic SEIR module for human infection to include three additional recovered states, each representing a boost in immunity (which is then lost to lower recovered echelons and then to the base susceptible compartment via first-order kinetics), based on Niger and Gumel (2008). Additionally, the model considered seasonally varying temperature profiles for different regions of sub-Saharan Africa. They predicted malaria burden (as measured in terms of the total number of new cases of infection) to increase with temperature in the range 16–28 $^{\circ }\hbox {C}$, but to decrease for temperature values above 28$^{{\circ }}$C in West Africa, 27 $^{{\circ }}$C in Central Africa, 26 $^{{\circ }}$C in East Africa and 25 $^{{\circ }}$C in South Africa. They also found that omitting either immature mosquito dynamics or temperature variability could significantly affect predictions, but the immunity-boosting module had little effect on infection incidence.

A similar effort by Okuneye and Gumel (2017) included partial resistance after infection, but also subdivided the human population into those under five (which bear, by far, the brunt of disease) and those over, and considered rainfall (adopted from Parham and Michael (2010)) in addition to temperature. This work found, for the Kwa-Zulu Natal province of South Africa, an increase in malaria burden with increasing mean monthly temperature and rainfall in the ranges 17–25 $^{\circ }\hbox {C}$ and 32–110 mm, respectively, and that malaria transmission is maximized for mean monthly temperature and rainfall in the ranges 21–25 $^{\circ }\hbox {C}$ and 95–125 mm. This model demonstrated dynamics only marginally affected by the inclusion of immunity and age-structure, in terms of infection incidence. While the notion that infection per se is minimally affected by age-structure and immunity is reasonably congruent with epidemiologic data, e.g. Trape et al. (1994), elucidating the effect of climate on serious clinical disease incidence (i.e., that disease which is severe enough to present clinically, as opposed to simply the acquisition of new infection, which may be asymptomatic among immunes) should be a future modeling goal. It should also be noted that, as in Mordecai et al. (2013), these works treat some interdependent processes as independent, e.g. the temperature-dependent biting and oviposition rates.

A parallel series of works by Yamana, Bomblies, and colleagues (Bomblies et al. 2008; Yamana et al. 2013, 2016), founded upon the agent-based HYDREMATS model developed by Bomblies et al. (2008) incorporated detailed hydrodynamic modeling at the village scale, and in a 2013 work (Yamana et al. 2013), the model framework was extended to include the gradual acquisition of partial immunity in humans, whereby repeated infections both reduced the probability of infection and increased the infection clearing rate, and concluded that, via immunity, large differences in infectious biting rate did not, in two highly malarious villages, translate into comparable differences in infection. This framework was also applied in a later (Yamana et al. 2016) model-driven effort that concluded climate change is unlikely to appreciably affect malaria burden in Western Africa.

6.3.5 Christiansen-Jucht model (2015) and age-dependent survival

It has frequently been observed that neither larval nor adult survival times are exponentially distributed (see Sects. 5.2.3 and 5.2.2), and this fact was incorporated into an age-structured vector lifecycle model by Lunde et al. (2013b), but the formulation was rather complex; note that delay-differential models and those with multiple age compartments (e.g., for larvae) also yield at least some non-exponential survival times.

We focus here on a more recent and easily digestible effort by Christiansen-Jucht et al. (2015) that included age-dependent survival at both the larval and adult stages (in addition to temperature-dependent survival), with data and model linked via the gamma distribution. The model entails an ODE framework with larval and adult populations divided into multiple age categories, with first-order transitions through categories and the final category terminating in death, as depicted in Fig. 17. We subdivide into four model combinations: (1) Baseline model without age-dependent survival, (2) Baseline + Larval age-dependent survival, (3) Baseline + Adult age-dependent survival, and (4) Larval and adult age-dependence.

Fig. 17

Generic technique for modeling age-dependent events in any population (left) and 2015 Christiansen-Jucht et al. model framework (right) employing this technique for the A. gambiae lifecycle, with age-dependent survival in both larval and adult populations

Fortuitously, mosquito (and other age-dependent) survival data can be reasonably well-described by the two-parameter gamma distribution, with $\alpha$ the shape parameter and $\beta$ the rate parameter, and it so happens that an ODE with $\alpha$ age compartments that are traversed at rate $\beta$ yields an overall gamma-distributed survival time (Wearing et al. 2005) (obviously $\alpha$ must be an integer for this to hold); this basic fact also facilitates ODE modeling of other non-Poisson time-dependent processes, such as infection clearance (Wearing et al. 2005). Note that if $\alpha$ = 1, then we reduce to exponentially distributed survival times. If $y_i(t)$ is the number of mosquitoes in compartment i, we have simply

$$\begin{aligned} \frac{dy_1}{dt}= & {} -\beta y_1, \end{aligned}$$

(73)

$$\begin{aligned} \frac{dy_i}{dt}= & {} \beta y_{i-1} - \beta y_i, \end{aligned}$$

(74)

and the total number of surviving mosquitoes at any time is

$$\begin{aligned} \sum _{i=1}^{\alpha } y_i(t). \end{aligned}$$

(75)

An example of a single cohort of adult mosquitoes is given in Fig. 18, using $\alpha$ = 7 and $\beta$ as fitted to adult survival data at 20 and 30 $^{\circ }\hbox {C}$ extracted from Bayoh (2001).

Fig. 18

Age-dependent death modeled using a multi-compartment model with $\alpha$ compartments and transition rate $\beta$, yielding a Gamma ($\alpha$,$\beta$) distributed survival curve. A seven-compartment model, i.e. $\alpha$ = 7, with $\beta$ fit to adult A. gambiae survival data from Bayoh (2001) at different temperatures. The left panels show how the distribution among age-compartments shifts over time, co-plotted with overall model survival and data, using $\beta$ for 20 and 30 $^{\circ }\hbox {C}$. The right panel gives $\beta$ as a function of temperature across the full data’s full temperature range

Now, using the framework given in Fig. 17, $\alpha$ and temperature-dependent $\beta$ parameters were fit for larvae and adults ($\alpha$ = 7 for larvae, $\alpha$ = 3 for adults), other temperature-dependent survival and larval development rates were adapted from Parham et al. (2012), and larval density-dependent death with carrying capacity a function of rainfall as adapted from White et al. (2011), Christiansen-Jucht et al. (2015) compared the ability of the four model permutations to fit mosquito abundance data from The Gambia, finding that Models 3 and 4 (adult age-dependent survival with and without larval age-dependent survival, respectively) gave nearly identical fits that were superior to either Model 1 or 2 (baseline or larval age-dependence only), suggesting that non-exponential death rates in adult mosquitoes are important for model fidelity to data.

6.3.6 Temperature variability and extreme weather

While most modeling works have considered ambient temperature as a single constant, in the last few years there has been increasing experimental and theoretical interest in the role of temperature variability, especially diurnal variation, on vector survival and development and malaria transmission potential (Paaijmans et al. 2009, 2010; Gething et al. 2011; Paaijmans et al. 2013b; Blanford et al. 2013; Lyons et al. 2013; Murdock et al. 2016; Beck-Johnson et al. 2017). For any given mean temperature point, survival and development rates do not generally change symmetrically with temperature excursions in either direction (given the non-linear nature of the related thermal-response functions), and thus we may expect exposure to fluctuating temperatures (as experienced in the field) to have a fundamentally different effect on vectors and parasites than exposure to a constant temperature.

This was first addressed in a theoretical context by Paaijmans et al. (2009), who showed that, using a unimodal Briere function (see Eq. (72)) for the Plasmodium sporogonic duration, fluctuations about relatively low temperatures enhanced development and hence malaria potential, relative to a constant temperature, while fluctuations at higher temperatures had the opposite effect, suggesting that most existing theoretical works may have systematically under- and overestimated malaria potential at cool and high temperatures, respectively. This pattern was subsequently observed experimentally in a rodent malaria model (Paaijmans et al. 2010), and Blanford et al. (2013) scaled these results up geographically to Kenya, predicting the cooler highlands to be relatively more vulnerable to malaria once diurnal temperature variation is accounted for; similar predictions were made at a continental scale as well. Paaijmans et al. (2013b) also found temperature fluctuations to enhance and inhibit A. stephensi development and survival at low and high temperatures, respectively, and similar results have been obtained for other ectotherms, e.g. Bayu et al. (2017).

The aforementioned works focused either on the sporgonic cycle (Paaijmans et al. 2009, 2010; Blanford et al. 2013) or Anopheles development (Paaijmans et al. 2013b) in isolation, but very recently, Beck-Johnson et al. (2017) explored daily and annual temperature variations under their previously discussed 2013 model (Beck-Johnson et al. 2013) (Sect. 6.3.1). This work suggests more subtle effects of temperature variation on both total adult mosquito and potentially infectious adult mosquito populations. Of particular interest, increasing the daily temperature range narrowed the mean temperature range over which both such populations could exist, by decreasing mosquito abundance at both higher and lower temperatures. That is, while the works above, which focused on a single aspect of the malaria transmission cycle suggested that temperature fluctuations asymmetrically favor transmission at low temperatures, the more complete model of Beck-Johnson et al. (2017) contradicts this notion to some degree. However, further complicating the picture, greater annual temperature variation tended to increase both the range of the infectious mosquito population (i.e. there were more times with both fewer and greater numbers of mosquitoes) and the mean at lower temperatures, while also decreasing the infectious population at the high temperature range. Thus, temperature fluctuations at different scales may affect malaria transmission in a variety of subtle ways that are not predictable from isolated components of the malaria lifecycle.

Several other models have incorporated seasonal temperature variability, e.g. Agusto et al. (2015), although this was not their primary focus. Also of note, Gething et al. (2011) using a Ross–Macdonald-style expression for vector potential, developed global maps for malaria suitability by temperature that was novel in that it superimposed diurnal sinusoidal temperature variations onto monthly temperature trends drawn from the WorldClim database, while a similar effort by Garske et al. (2013) inferred air temperatures from land temperature data and also imposed diurnal temperature variation on a Ross–Macdonald-like formulation for malaria potential. We suggest that comparing the predicted malaria maps using more complex models with and without such temperature variations would be a valuable contribution to the literature.

It should be mentioned that the works reviewed here have generally considered ambient air temperature only. Yet, as already discussed, the temperature in aquatic anopheline habitats may differ appreciably from the air, and, depending upon the habitat size and thermal stability, temperature variations in water may be smaller or greater than in air, and how this complexity might further alter the predicted role of temperature fluctuations, both current and as anticipated under climate change, is an open question. Notably, the work by Blanford et al. (2013) did consider indoor temperatures, which tend to be slightly higher but less variable than ambient (Singh et al. 2016; Blanford et al. 2013), and the indoor/outdoor temperature difference has been almost uniformly neglected in modeling works (but see Singh et al. (2016) for a recent exception), despite the indoor preference seen in many anophelines.

Finally, climate change is likely to increase weather extremes, including drought, extreme rainfall events, and heat waves, which may be expected to affect the lifecycles of a range of ectotherms, independent of average weather conditions (Ma et al. 2015). Southern Africa, especially, may be vulnerable to greater frequency of extreme rainfall events (Engelbrecht et al. 2013). While a potential increase in immature Anopheles mortality due to heavy rainfall is sometimes accounted for (Paaijmans et al. 2007) to our knowledge no mathematical work has more explicitly examined how weather extremes under climate change, and especially extreme high temperatures, might affect malaria epidemiology.

6.4 Towards a meta-population model for the Kenyan highlands

Recall that the Kenyan highlands have seen malaria incidence increase since the 1970s in conjunction with increasing temperatures and broad changes in the populace, including rapid population growth and deforestation, making this a model region for studying the impact of current and projected climate change on malaria transmission (Minakawa et al. 1999; Pascual et al. 2006; Chaves and Koenraadt 2010). Human mobility can link regions with varying local transmission dynamics (e.g. as a result of climate), and while no climate-focused models have accounted for this thus far, we present basic mathematical frameworks by which this might be studied in the future. Any geography may be conceptually divided into multiple patches, each with its own sub-model describing local disease dynamics, and with movement of hosts and/or vectors between patches occurring at prescribed rates; thus movement of Plasmodium reservoirs from one patch to another may spread disease. A natural division is to consider highland and lowland areas as two separate patches, but any two-patch model easily generalizes to an n-patch model.

It is necessary to distinguish between two forms of mobility: (1) migration, representing permanent resettling in a new area, and (2) visitation, or transient excursions with no permanent change of address (Mandal et al. 2011). Both are salient, as large-scale migration into the Kenyan highlands has occurred over the past few decades (Minakawa et al. 1999), while small-scale and circulatory movements between communities are also common. Indeed, transient mobility also interacts with age because, in many parts of rural Africa, mothers from rural areas tend to take their infants and young children (not of school-going age) on a short trip (usually for a day or two) to conduct businesses at open markets in neighbouring urban communities, thereby exposing them to malaria infection, especially if this trip is from highland areas of low endemicity to lowland areas of high malaria burden. In a multi-patch model, permanent migration is represented, in analogy to fluid dynamics, by the classical Eulerian approach, whereas transient mobility with a home patch is captured by the Lagrangian approach (Castillo-Chavez et al. 2016).

Mathematically, Eulerian migration resembles a diffusion process, and we define $k_{ji}$ to be the first-order rate-constant for movement from patch j to i. A very simple Ross-style model with Eulerian migration among human hosts, with $H_i$ the total human population and $X_i$ the infected population in patch i (and $M_i$ and $Z_i$ the total and infected mosquito populations), is

$$\begin{aligned} \frac{dH_i}{dt}= & {} \sum _{j=1, j \ne i}^\varPhi k_{ji} H_j - \sum _{j=1, j \ne i}^\varPhi k_{ij} H_i, \end{aligned}$$

(76)

$$\begin{aligned} \frac{dX_i}{dt}= & {} a b \left( \frac{Z_i}{H_i}\right) (H_i - X_i) - r X_i + \sum _{j=1, j \ne i}^\varPhi k_{ji} X_j - \sum _{j=1, j \ne i}^\varPhi k_{ij} X_i, \end{aligned}$$

(77)

$$\begin{aligned} \frac{dZ_i}{dt}= & {} a c \left( \frac{X_i}{H_i}\right) (M_i - Z_i) - g Z_i, \end{aligned}$$

(78)

where $\varPhi$ is the total number of patches. In other words, it is simply the standard model augmented by migration terms among the hosts. The Lagrangian formulation is slightly less obvious. We have $p_{ij}$ as the fraction of time that individuals from patch i spend in patch j, subject to the constraint that $\displaystyle \sum \nolimits _{j=1}^\varPhi p_{ij} = 1$. We also have that the effective population of patch i is $\displaystyle \sum \nolimits _{k=1}^\varPhi p_{ki} H_k$. A Ross model with Lagrangian motility, between an arbitrary number of patches, is therefore given by

$$\begin{aligned} \frac{dX_i}{dt}= & {} a b \displaystyle \sum _{j=1}^\varPhi \left( {\frac{p_{ij} Z_j}{{\displaystyle \sum _{k=1}^\varPhi (p_{kj} H_k)}}}\right) (H_i - X_i) - r X_i, \end{aligned}$$

(79)

$$\begin{aligned} \frac{dZ_i}{dt}= & {} a c \displaystyle \sum _{j=1}^\varPhi \left( \frac{p_{ji} X_j}{p_{ji} H_j}\right) (M_i - Z_i) - g Z_i. \end{aligned}$$

(80)

Note that hosts do not move between home patches, but the effective population of each patch is modulated by the time fraction every other patch population spends within it. The two mobility modes can easily be combined into a single model, essentially by augmenting Eqs. (79)–(80) with the Eulerian mobility terms.

Multiple authors have studied multi-patch Ross–Macdonald style models (Torres-Sorando and Rodríguez 1997; Auger et al. 2008; Cosner et al. 2009; Prosper et al. 2012; Ruktanonchai et al. 2016), beginning with Eulerian mobility in Torres-Sorando and Rodríguez (1997) and Lagrangan mobility in Cosner et al. (2009); see also Agusto (2014) for a brief review. Agusto (2014) recently studied a more complex (but weather-independent) multi-patch malaria model that also included the vector lifecycle, and moreover, focused on the spread of drug-sensitive and drug-resistance Plasmodium strains under Eulerian migration; to our knowledge, mobility has not been incorporated into a weather-driven model using the patch framework, and we suggest it as an important extension of current models. Additionally, we offer our speculation that a multi-patch model at very local scales with mobility among vectors, in addition to (or in lieu of) hosts, could help elucidate how the presence of varied microenvironments across a fine spatial scale might alter transmission dynamics.

7 Other modeling challenges: immunity, treatment, within-host disease, and other abiotic factors

In this section, we briefly touch on some other malaria modeling challenges and traditions, including immunity and within-host dynamics, treatment, resistance, and socioeconomic factors. While there is insufficient space do any of these topics justice, we believe that incorporating many of these aspects of disease in climate-driven model frameworks is a fundamental challenge for the future, and hope to at least make the reader aware of these issues and some useful references.

7.1 Malaria immunity

Numerous studies have attempted to quantify the acquisition of protection against clinical malaria with repeated infection, beginning as early as the Garki model of Dietz et al. (1974) (Sect. 4.4), and include (Aron 1983, 1988; Gupta and Day 1994; Gupta et al. 1999a, b; Filipe et al. 2007; Griffin et al. 2010, 2015). Filipe et al. (2007), for example, concluded that a short-term, primarily clinical, immunity strongly coupled to cumulative exposures and a longer-term, primarily anti-parasitic, immunity more weakly coupled to exposure are both necessary to best fit epidemiologic data. These authors modeled infection using an age-structured model, where exposed persons either manifest clinical disease or asymptomatic disease, with the probability of clinical disease modulated by an underlying level of immunity, $I_s$, which increases in direct proportion to EIR, and decays exponentially with a decay half-life of about 5 years. The probability of clinical (versus asymptomatic) disease, $\phi$, decreases according to a sigmoid function of $I_s$. Additionally, asymptomatic infection is assumed to clear at an exponential rate that increases with age, given that some exposure has occurred, but otherwise independent of cumulative exposure; the half-life of this parasitic immunity was determined to be at least 20 years. A similar approach was also used recently by Griffin et al. (2010) and Griffin et al. (2015).

This model-based conclusions is concordance with clinical data per Rodriguez-Barraquer et al. (2016), who performed an analysis of detailed longitudinal data of 93 children over the first 5 years of life in a holoendemic region of Uganda, who underwent surveillance for both clinical disease and (microscopy-detected) asymptomatic parasitemia. Infection was modeled essentially as a Poisson process, with individual infection hazard varying according to a Gamma distribution.

This study suggested a binary effect of age, with increasing age a risk factor for infection, but also independently protective against malaria, given infection (regardless of infection history). The models suggested a 6% decrease in malaria given infection per year of age, and a 2% decrease in malaria per infection. Note, however, that there were 5.33 episodes of malaria per person-year (and 0.588 asymptomatic parasitemias per person-year), suggesting that infectious history was dominant over age in conferring protection. A further finding was that clinical immunity developed faster in children treated with artemether-lumefantrine (AL) versus dihydroartemisinin-piperaquine (DP). The latter drug confers longer-lasting post-treatment protection against infection.

As already alluded to multiple times, given its apparent historical and epidemiologic importance, further mathematical investigations incorporating immunity with climate factors are essential. Towards that end, several recent climate-focused works have included relatively simple representations of partial immunity (Agusto et al. 2015; Yamana et al. 2013, 2016, 2017), although they have not generally considered the disparity between anti-disease and anti-parasite immunity. In particular, work by Yamana et al. (2013, 2017) represents important steps towards understanding the interplay between environment and immunity.

7.2 Treatment, control measures, and resistance

There was a large drop in the African malaria burden in the 1960s and 70s, largely attributable to widespread pharmacologic treatment with chloroquine, but this was undermined by the evolution and spread of chloroquine resistance P. falciparum strains (Carter and Mendis 2002). Artemisinin compounds are now highly effective, but resistance has already been detected in southeast Asia, and it seems that the basic biology dictates that it is only a matter of time before resistance reaches Africa (Webb 2014). A wide variety of mathematical models have addressed the evolution of drug resistance in both infectious disease and cancer. Some works focused on modeling drug resistance in malaria include (Hastings 1997, 2003; Yeung et al. 2004; Pongtavornpinyo et al. 2009; Saralamba et al. 2011; Agusto 2014; Forouzannia and Gumel 2015), and Hastings and Watkins (2005) provide an excellent review of the main concepts involved in modeling resistance. These are only a few examples from the literature, but any deeper review of these works is unfortunately beyond our scope.

We also note that numerous models have focused on other control efforts, such as ITNs, (most classically we have Macdonald’s prediction that adult mosquito survival should be targeted, motivating insecticide-based control), with some recent efforts including Smith et al. (2009), Griffin et al. (2010) and Nikolov et al. (2016), and finally, a recent and broad review of the evolutionary principles underlying resistance is given by Huijben and Paaijmans (2017).

7.3 Within-host disease dynamics

Coupling an explicit model for the within-host dynamics of malaria (principally the within-host immune response) to its epidemiology is a fundamental challenge. Indeed, even describing the within-host dynamics mathematically in a way the reproduces the qualitative dynamics of long-term infection has proven most difficult, and this issue has its own extensive literature that it is beyond the scope of this paper, although a very partial reference list includes Teboh-Ewungkem et al. (2010), Li et al. (2011), Saralamba et al. (2011), Gurarie et al. (2012), Eckhoff (2012), Demasse and Ducrot (2013), Childs and Buckee (2015), Childs and Prosper (2017), Tabo et al. (2017), and a recent work by Childs and Buckee (2015) highlights the problems of accurately modeling within-host dynamics in some depth. Since these dynamics may interact in unexpected ways with malaria epidemiology (Childs and Buckee 2015), a complete understanding of climate and malaria will likely necessitate a deeper consideration of the in-host stages than has heretofore been attempted.

7.4 Other abiotic factors

Broad changes in socioeconomic conditions, e.g. widespread urbanization and increasing material standards of living, and agricultural modernization played fundamental roles in the retreat of malaria from most of the world outside tropical Africa (Webb 2014; Packard 2007). Socioeconomic conditions also determine access to medical treatment and effective antimalarials (Webb 2014), and rural populations suffer a higher malarial burden than do urban populations (Rodriguez-Barraquer et al. 2016), partly because urban land-use patterns support fewer vectors and partly because control efforts have historically focused on urban over rural areas (Webb 2014). Land-use changes, driven by social or economic imperatives, strongly affect anopheline habitat, e.g. deforestation in the Kenyan highlands (Afrane et al. 2005). Incorporating social factors directly into models is a challenge, but can probably be expected to affect biophysical parameters in somewhat predictable downstream ways; see also Mandal et al. (2011) for a brief review of models relating socioeconomics to malaria.

8 Conclusions

Of the historical human Plasmodia, P. falciparum is evolutionarily distant from the others (Silva et al. 2015), and uniquely virulent. While in the West, the so-called diseases of civilization are generally thought of as those cardiovascular ailments, diabetes, and obesity that tend to accompany high-energy diets and sedentary lifestyles, P. falciparum malaria was perhaps the first true disease of civilization, ushered in by global warming, anthropogenic alteration of the environment, and concentrated human settlement. Malaria has been subject, more than most diseases, to mathematical investigations, and these efforts, mainly the pioneering works of Ross and Macdonald, had real influence on malaria control efforts through the twentieth century. It is likely, then, that the mathematician can play a central role in informing an understanding and mitigation of global warming’s future impact of malaria. Malaria, however, is a complex disease with its own particular history, and so we take the view that mathematical efforts are best informed by history, and to that end have attempted a reasonably thorough historical review of the disease.

Globally, malaria retreated dramatically over the course of the twentieth century, most markedly in post-World War II Southeast Asia (Carter and Mendis 2002). This drop occurred in the face of very modest global warming (about 0.6 $^{\circ }\hbox {C}$) (IPCC 2013; Gething et al. 2010), and is almost certainly attributable to a variety of non-climatic changes, especially economic and agricultural modernization, urbanization, and broad increases in population health and health services (Chaves and Koenraadt 2010; Webb 2014). While this would seem to suggest climate change as a minor, at best, factor in future transmission scenarios, the consensus view of the IPCC is that risks from climate change are not likely to be strongly felt until the global temperature anomaly exceeds at least 1 $^{\circ }\hbox {C}$, with impacts increasing dramatically above 2 $^{\circ }\hbox {C}$ (IPCC 2014); the temperature anomaly is virtually guaranteed to reach 1.5 $^{\circ }\hbox {C}$ by the end of the twenty-first century, and may well exceed 5 $^{\circ }\hbox {C}$ (IPCC 2014).

Tropical Africa never saw the deep reductions in malaria mortality that other regions did and, while it enjoyed a transient drop in deaths in the 1960s and 70s, malaria re-surged in the later 1970s, an era corresponding with widespread chloroquine resistance, other broad socioeconomic changes (Carter and Mendis 2002; Webb 2014), and the acceleration of global warming (IPCC 2014). While recent control efforts have been reasonably successful, with malaria mortality in Africa now at an historical nadir, their sustainability is threatened on multiple fronts (see Sect. 3.3), and future temperature increases in Africa are likely to exceed the global mean (Niang et al. 2014), where global warming may negatively affect agricultural output, food security, economic development and overall health, and may displace populations (Niang et al. 2014), all developments likely to increase populations’ vulnerability to malaria. Thus, it seems likely that future warming will interact with malaria (and other infectious disease) in a nonlinear manner. Therefore, we must be cautious in extrapolating from past climate and malaria trends to the future, and a mechanistic framework for the disease and climate, based on the many excellent works reviewed in this paper, may help to guide us.

There are perhaps five major challenges to mechanistically modeling the relationship between climate change and malaria: (1) thermal-response functions linking temperature and vector/lifecycle parameters; (2) the relationship between weather (mainly rainfall) and anopheline habitat availability, and further, the effect of habitat size on the vector; (3) temperature variability, at both a diurnal and seasonal scale, and in the different microenvironments involved in the Anopheles lifecycle; (4) the incorporation of essential non-climatic factors, especially malaria immunity, but also treatment and other control interventions, host mobility across zones of varying endemicities, resistance, and broad socioeconomic factors; and finally (5) basic model construction, i.e. the general choice of biologic actors and their interactions. Most controversy has more explicitly centered on the first two, especially thermal-response functions, and to that end we have presented a detailed summary of experimental data to inform the modeler (Sect. 5), as well as a partial genealogy of the work informing this controversy. Box 1 provides a summary of some of the most important, at least in our view, modeling challenges moving forward.

While the earlier works of Martens and colleagues (e.g. Martens et al. (1999)) suggested a significant global increase in potential malaria burden with global warming, these efforts were informed by thermal-response functions that fail to capture the deleterious effects of high temperatures on both vector and parasite. The unimodal thermal-response functions of Mordecai et al. (2013) suggest a lower temperature range for optimal transmission (25–28 vs. 32–33 $^{\circ }\hbox {C}$), but also view some interconnected components of the vector lifecycle as independent, e.g. biting rate and oviposition, and death prior to larval development may be conflated with arrested development. Mechanistic models have properly captured these lifecycle interdependencies to lesser and greater degrees (e.g. Hoshen and Morse 2004; Parham and Michael 2010), and while several later works employing the Mordecai et al. relations in more complex mechanistic models (Agusto et al. 2015; Okuneye and Gumel 2017), have reached conclusions generally congruent with Mordecai et al. (2013), these too view certain dependent processes as independent. Therefore, more careful inclusion of thermal-response functions in mechanistic frameworks should be a goal of future work.

Relating rainfall to the vector lifecycle is less prominent in this modeling tradition, and when it has not been ignored entirely, its presumed effect is more variable and ad hoc across models. Nevertheless, precipitation patterns are likely to change across Africa with climate change, and rainfall patterns often drive interannual malaria variability (Pascual et al. 2008). More realistic hydrodynamic models are likely to be informative, and we suggest the hydrodynamic models of, for example, Bomblies et al. (2008), Parham et al. (2012) and Asare et al. (2016a) as a starting point; some of these models’ basic features and construction are reviewed in Sect. 5.2.5.

Basic model construction, including weather-independent components, is clearly fundamental, with published models encompassing a wide range of biological detail and realism. Mechanistic works have demonstrated that greater mathematical fidelity to the details of the vector lifecycle, e.g. via the inclusion of immature larval stages (Agusto et al. 2015; Beck-Johnson et al. 2013) or age-dependent survival (Christiansen-Jucht et al. 2015), significantly improves model predictions in relation to data. Thus, it is important to elucidate how such deeper model construction choices affect predictions, in conjunction with various thermal- and rainfall-response functions. Furthermore, malaria immunity is fundamental to its epidemiology, and while largely neglected in climate-focused works, there is a large body of immune-focused models, and several recent works have demonstrated that immunity interacts importantly with environment (Yamana et al. 2013, 2017).

Human mobility, both via long-term migration and shorter-term circulation of individuals between areas of low and high malaria burden, may affect malaria incidence, but it is unknown as yet how this mobility dynamic interacts with climate, and we suggest that a multi-patch model that takes into account the detailed lifecycles of vector and parasite, temperature variability and the physics of anopheline habitat, along with the host phenomena of immunity, superinfection, and asymptomatic infection would entail a somewhat unified framework for studying malaria spread from lowland to highland regions. We further suggest that this more limited geographic scope may better elucidate the effect of climate change on malaria transmission than global scale models. Nevertheless, the development of a new family of malaria potential maps, employing a variety of more recent malaria models, temperature variability, and IPCC climate projections, may be highly instructive in determining the robustness of past authors’ conclusions and sensitivity to modeling choices. Towards such an end, the recent works of Ryan et al. (2015b) and Yamana et al. (2016) can serve as excellent guides.

Additionally, making the primary focus of malaria mapping upon populations at risk, both current and projected, rather than land area, is likely of import, especially since most population in Africa is currently concentrated in relatively warm Western coastal Africa, where malaria is highly endemic and the effect of climate change upon malaria may be equivocal, and in cooler areas of Eastern Africa, mainly Ethiopia and the region surrounding Lake Victoria, where global warming may be more likely to increase disease potential. Dramatic population growth is projected throughout sub-Saharan Africa, and Nigeria, already the most populous African nation and one with a high malaria burden, may see its population more than double by 2050 (United Nations 2017). Urban environments also are less malarious than the countryside, and so the urban versus rural divide is also important to consider (as done in Ryan et al. (2015b)).

There are multiple modeling traditions within the field of mathematical malariology, and we have been restricted by space to focus primarily on only one; we have very briefly touched on these traditions in Sect. 7, and we also refer the reader to other useful reviews, including Mandal et al. (2011), Smith et al. (2012) and Reiner et al. (2013). Malaria immunity, abiotic factors including land use, and control efforts have all proven fundamental to the epidemiology of the disease, and unifying weather-driven models with these biologic and social phenomena is a fundamental task for the future.