- Research
- Open Access
- Published:

# Using Bayesian state-space models to understand the population dynamics of the dominant malaria vector, *Anopheles funestus* in rural Tanzania

*Malaria Journal*
**volume 21**, Article number: 161 (2022)

## Abstract

### Background

It is often assumed that the population dynamics of the malaria vector *Anopheles funestus*, its role in malaria transmission and the way it responds to interventions are similar to the more elaborately characterized *Anopheles gambiae*. However, *An. funestus* has several unique ecological features that could generate distinct transmission dynamics and responsiveness to interventions. The objectives of this work were to develop a model which will: (1) reconstruct the population dynamics, survival, and fecundity of wild *An. funestus* populations in southern Tanzania, (2) quantify impacts of density dependence on the dynamics, and (3) assess seasonal fluctuations in *An. funestus* demography. Through quantifying the population dynamics of *An. funestus*, this model will enable analysis of how their stability and response to interventions may differ from that of *An. gambiae *sensu lato.

### Methods

A Bayesian State Space Model (SSM) based on mosquito life history was fit to time series data on the abundance of female *An. funestus *sensu stricto collected over 2 years in southern Tanzania. Prior values of fitness and demography were incorporated from empirical data on larval development, adult survival and fecundity from laboratory-reared first generation progeny of wild caught *An. funestus*. The model was structured to allow larval and adult fitness traits to vary seasonally in response to environmental covariates (i.e. temperature and rainfall), and for density dependency in larvae. The effects of density dependence and seasonality were measured through counterfactual examination of model fit with or without these covariates.

### Results

The model accurately reconstructed the seasonal population dynamics of *An. funestus* and generated biologically-plausible values of their survival larval, development and fecundity in the wild. This model suggests that *An. funestus* survival and fecundity annual pattern was highly variable across the year, but did not show consistent seasonal trends either rainfall or temperature*.* While the model fit was somewhat improved by inclusion of density dependence, this was a relatively minor effect and suggests that this process is not as important for *An. funestus* as it is for *An. gambiae* populations.

### Conclusion

The model's ability to accurately reconstruct the dynamics and demography of *An. funestus* could potentially be useful in simulating the response of these populations to vector control techniques deployed separately or in combination. The observed and simulated dynamics also suggests that *An. funestus* could be playing a role in year-round malaria transmission, with any apparent seasonality attributed to other vector species.

## Background

*Anopheles funestus* is one of the major malaria vectors in Africa and is widely distributed across the continent [1, 2]. With the exception of *Anopheles gambiae *sensu stricto (*s.s*.), the species appears to have higher vectorial capacity than many other members of the *Anopheles gambiae* complex [3,4,5,6,7,8]. *Anopheles funestus* makes a higher contribution to transmission than *An. gambiae *sensu lato (*s.l*.) in numerous parts in sub-Saharan Africa [6, 9,10,11,12]; particularly in settings where *An. gambiae* abundance has plummeted due to either effective indoor-based vector control interventions [13, 14] or environmental change. It is hypothesized that *An. funestus* persistence despite the recent scale-up of insecticide–treated nets may have been facilitated by their earlier development of strong physiological resistance [15].

*Anopheles funestus* is typically grouped with *An. gambiae s.l.* when modelling transmission and formulating policies for malaria vector control [16,17,18]. The lack of explicit consideration of *An. funestus* ecology and transmission potential may be partially due to this species having been relatively neglected compared to *An. gambiae*. Comparatively the ecology of *An. funestus s.l*. is less well understood, and it is much more difficult to maintain under insectary or semi-field conditions [19]. However, this species has several unique ecological features, such as its different larval habitat and dry season persistence [20], that could give rise to distinct population dynamics and differentiate its response to core and supplementary interventions. For example, *An. funestus* prefers larger aquatic habitats that are semi-permanent or permanent throughout the year, and contain clear water with some emergent vegetation [20]. This differs from *An. gambiae s.s*. which generally prefer small temporary habitats, such as puddles, ditches or animal hoof prints [20, 21], or *Anopheles arabiensis*, which can breed extensively in rice fields and other sunlit open pools [1]. The use of more permanent larval habitats means that *An. funestus* has greater persistence through the driest periods of the year compared to *An. gambiae* [22], whose habitats evaporate quickly in the absence of rainfall [21, 23]. This ecological feature means that the seasonal phenology of *An. funestus* and its response to aquatic microclimate differs from *An. gambiae* [21, 22, 24]; and could thus generate differential response to seasonally-targeted interventions, such as Indoor Residual Spraying (IRS) and larviciding.

Differential use of aquatic habitats may also impact the relative importance of key intrinsic drivers of mosquito population dynamics such as density dependence. Density dependence in malaria vectors occurs during larval development as a product of competition for space and nutritional resources [25, 26]. In space-limited habitats, high larval densities can influence larval development rates and survival, but also subsequent adult fitness traits such as body size, survival, fecundity and mating success [27,28,29,30]. While there is evidence that density dependence is an important driver of *An. gambiae* population dynamics [25], the relative importance of this process for *An. funestus* is unknown. Given that larval crowding and competition are less likely within the larger habitats preferred by *An. funestus,* density dependence is hypothesised to may be less pronounced for this vector species. Quantifying the strength of density dependence is important to inform the ease with which vector populations can be suppressed and how quickly they can recover [26, 27, 31].

Models of vector population dynamics and their response to interventions must be parameterised by reliable estimates of their demography and fitness. For vectors in the *An. gambiae* complex, such estimates are often acquired from insectary and semi-field studies [32,33,34,35,36,37] as well as field studies. Similar data has been difficult to obtain for *An. funestus* because of its poorly understood ecology and the difficulties of creating laboratory colonies; which so far has been achieved on only two occasions [11, 19, 38]. State-space models (SSM) provide an alternative approach to indirectly estimate these parameters by fitting a population dynamics model to observed time series data [39, 40]. These models are widely used in other fields of ecology and conservation biology to investigate the population dynamics of other animals [39, 41] and guide management decisions [41]. However, these models have so far had limited update in medical entomology. Given data on population fluctuations are available, these models can infer and estimate plausible demographic rates that could generate the observed dynamics [42].

SSMs are time-series models that distinguish between two stochastic components, namely, process (i.e. biological), which captures sequential dependencies between population components (e.g. eggs, larvae, pupae) and an observation component, which captures and corrects for biases and imprecisions in the data-collection process. Prior knowledge of the model parameters is used to bolster the information content of the time series data with existing expert or laboratory data and uncertainty in estimates. Population projections are then quantified on the basis of posterior probability distributions for parameters and population states. SSMs have recently been used to elucidate the dynamics and impacts of interventions on malaria vectors in laboratory and semi field populations [32, 33], but have not yet been applied to estimate *An. funestus* vector demographics in the wild. Here, an innovative SSM application was developed to describe the dynamics of wild *An. funestus* populations in Tanzania, and use it to assess extrinsic (environmental) and intrinsic (density dependence) drivers of their fitness and demography for the first time. Empirical data from laboratory experiments on *An. funestus* colonization [19] were incorporated together with the wild population data to develop an SSM. Time series field data collected in 2015 [6] and 2018 south-eastern Tanzania, and corresponding environmental information were used to validate the model. Specific aims were: (1) to accurately reconstruct the population dynamics, survival and fecundity of wild *An. funestus* populations in southern Tanzania, (2) quantify the effects of density dependence on the dynamics, and (3) to identify and quantify seasonal variations in *An. funestus* demography.

## Methods

### Time series data on wild *An. funestus* populations

Indoor densities of female *An. funestus s.l.* adults were recorded over 12 months of entomological surveys conducted in three villages (Tulizamoyo, Ikwambi and Sululu) in Kilombero (8.1539ºS, 36.6870ºE) and Ulanga (8.3124ºS, 36.6879ºE) districts, south-eastern Tanzania from June 2018 to May 2019 (Fig. 1). The villages were selected because of the high abundance of *An. funestus s.l.* within which *An. funestus s.s.* is the dominant sibling species (93%) [1]. Annual rainfall was 1200–1800 mm, and temperature, 20–32 °C. CDC Light traps [43] were used to sample host-seeking mosquitoes from 6 pm to 6am for 5 days per week, 4 weeks a month for 12 months in 10–15 houses per village. The houses were randomly selected and consent obtained from household heads, mosquito sampling were done repeatedly in these houses. The mosquitoes were sorted by taxa and sex, and females further classified as unfed, blood-fed or gravid. Daily climatic data (rainfall and temperature) were obtained from a weather station, approximately ~ 20 km from the farthest village.

To complement this, additional data on *An. funestus* were extracted from a 2015 dataset from three other villages in Ulanga district (Mavimba, Minepa and Kivukoni) (Fig. 1) [6]. These data were collected five days per week for a period of 12 months. This data allowed us to fit the model simultaneously to multiple time series so that it could learn hierarchically from *An. funestus* trajectories enfolding in different years and locations. This additional data has previously been described elsewhere and used to demonstrate the epidemiological dominance of *An. funestus*, which now contributes > 85% of all malaria infections in the region [6].

### Prior information on life-history and gonotrophic cycle stages

Female *An. funestus* adults collected from the same three villages in 2018 were maintained in insectary conditions for one generation to estimate baseline fitness traits as already described in Ngowo et al. [19] and in Fig. 2. Data collected from this 1^{st}-generation laboratory progeny included: (a) proportions of eggs that hatched into larvae, larvae that transitioned to pupae, and of pupae that emerged into adults, (b) the length of the transition periods (days) between life stages [(i) eggs to 1st instar larvae, (ii) 1st instar larvae to pupae, and (iii) pupae to unfed adult female (1 day post emergence)], (c) transition period of adult females between three different stages of their gonotrophic cycle, i.e. unfed, blood-fed and gravid.

The gonotrophic cycle starts with ‘unfed’ females who transition to ‘blood-fed’ after obtaining a blood meal. In the wild, the first gonotrophic cycle usually starts after unfed females have mated [44]; which is assumed to happen soon after emergence. In insectary experiments, females had access to males immediately on emergence. As the blood meal is digested, blood-fed females transition into the ‘gravid’ state during which eggs develop. Gravid females then oviposit their eggs into aquatic habitats and return to the ‘unfed’ stage with the cycle begins again (Fig. 2). In the wild, the rate of transition between these gonotrophic stages is governed by both intrinsic and extrinsic environmental conditions including the availability of blood-meals and oviposition sites [45]. In the insectary, the first blood-meal (arm-feeding) was offered 5 days post emergence to ensure individuals had sufficient time for mating.

Per capita fecundity was defined as the number of eggs laid per fully bloodfed adult female. The proportions surviving between life-history stages or gonotrophic stages were calculated as the inverse of the number of days required to transit from one stage to the next:

Here, \(m\) is the life cycle stage, \(j\) is proportion survived as the percentage of the total from preceding life stage, \(r\) is the average number of days it took to transit from one life stage to another, and \(s\) is the daily survival within the stage (Table 1).

### Biological process components of the Bayesian SSMs

#### Daily survival

The daily survival of larval stages was assumed to be the same for all instar stages. Adult survival was assumed to be the same for unfed, blood-fed and gravid females. Survival probabilities (\({s}_{l}, {s}_{a}, {s}_{f}, {s}_{v}\)) were linked to their covariates through a *logit* transformation of linear predictors (here, subscripts \(l, a, f, v\) refer to larvae, unfed, blood-fed and gravid, respectively). Pupal (\({S}_{p}\)) and egg (\({S}_{e}\)) survival probabilities were considered to be independent of any climatic and density-dependent covariates, and were treated as a binomial distribution, with baseline rates assigned priors as described in Eq. 16–17.

A range of covariates hypothesized to be associated with the demography of *An. funestus* were incorporated to allow baseline larval and adult survival to vary with environmental conditions. Rainfall (current and 1 week-lagged) and temperature were incorporated into the larval survival model. Rainfall regulates the availability and permanence of aquatic habitats, thus influencing both survival and carrying capacity of larval habitats [46]. Density dependence was incorporated into the model of larval survival [25] to assess whether this could improve the fit of the adult population dynamics model. Additionally, the speed of larval development was modelled as a function of temperature based on its known importance [47, 48]. The daily survival of larvae was thus defined as a function of daily rainfall (current and lagged), daily temperature and density dependence. The daily survival rates (lowercase \({s}_{l}\left(t\right)\)) of larvae were estimated through a logit transformation of linear predictors (uppercase \({S}_{l}\left(t\right)\)).

Specifically, \({S}_{l}\left(t\right)\) is written as a function of both intrinsic and extrinsic drivers:

Here, \({\beta }_{0}\) is the baseline daily larval survival on the linear scale. A survival probability prior was assigned under zero rainfall and average temperature (i.e. 27 °C) and then calculated the intercept of \({\beta }_{0}\) to reflect this prior information. When there is no effect of any environmental covariates (prior takes values between 0.80 and 1, Table 1). The coefficient \({\beta }_{1}\) quantifies the effect of current rainfall \((R)\)*;* with the envisioned scenario being that higher \(R\) (i.e. flooding) tends to wash away larvae hence reducing the baseline survival [49]. This \({\beta }_{1}\) was defined by an informative gamma prior with shape = 5.382 and rate of 46.4 (Table 2) which permits anything from no rain effects to 100% mortality. The coefficient \({\beta }_{2}\) quantifies the effect of larval density at time \(t\) on larval survival. A monotonic negative relationship was assumed based on the biologically-plausible hypothesis that larval survival is reduced at high larval density because of resource competition and intraspecific cannibalism [29, 50]. This coefficient \({\beta }_{2}\) was defined by an uninformative gamma prior with shape of 0.5 and rate of 1 (Table 2), which allows the impact of density to range from no effect to complete annihilation.

The term inside brackets in Eq. (3) represents the fact that density dependence needs to be modulated by the availability of larval habitat. The availability of suitable aquatic habitats for oviposition will increase with rainfall; thus potentially reducing the crowding of larvae into the remaining habitats that persist during the dry season. This hypothesis has been supported for *An. gambiae s.l*., where their seasonal population dynamics can be explained by models incorporating a rainfall-dependent carrying capacity [25]. Here, the coefficient \({\beta }_{3}\) was a proportion that captures the potential interaction between larval habitat availability (defined as the cumulative rainfall (Q) over the past week) and larval density (\(D\)). When rain in the recent week has been the maximum observed (i.e.\(Q=\mathrm{max}(Q))\), then (\(1 -\frac{{\beta }_{3}{Q}_{\left(t-1\right)}}{\mathrm{max}\left(Q\right)}\)) would be the smallest amount of density dependency experienced by *An. funestus*. The prior distribution for \({\beta }_{3}\) was defined by an upward-biased beta prior with mean 0.9 and variance of 0.01 allowing \({\beta }_{3}\) to have positive impact on larvae survival.

Additional covariates were incorporated to assess the role of temperature on larval survival (via the coefficients \({\beta }_{4}\),\({\beta }_{5}\)). The parameter \({\beta }_{4}\) captures the potentially positive effects of temperature on daily larval survival, which were defined by an uninformative gamma prior with mean of 1 and variance of 0.1 considering 27 °C as the optimal temperature (\(\rho )\) for maximum survival [51]. This prior allows temperature to vary from having no impact, to high positive impact on larval survival. Alternatively, the relationship between larval survival and temperature may be characterized by survival being reduced at low or very high temperature, and peaking in the middle [52]. The coefficient \({\beta }_{5}\) was incorporated to capture this potential curvilinear relationship but was dependent on \({\beta }_{4}\), to ensure that the optimum temperature was fixed at 27 °C, (\({\beta }_{5}=\frac{{\beta }_{4}}{2\rho })\). The prior of \({\beta }_{4}\) allowed extreme temperatures values away from the optimum to range from having no effect to generating 100% mortality. The parameter \({\varepsilon }_{l}\) is capturing the unexplained stochasticity associated with larval survival. This error term was defined by a normal prior with mean of 0 and a precision \({\sigma }_{l}\) from a gamma distribution with both shape and rate of 10.

The linear predictors for survival of unfed, blood-fed, gravid (uppercase \({S}_{a}\left(t\right), {S}_{f}\left(t\right), {S}_{v}(t)\)) and daily probabilities of survival (lowercase \({s}_{a}\left(t\right),{ s}_{f}\left(t\right), { s}_{v}\left(t\right)\)) were structured similarly to Eq. 3. The daily survival probabilities of adult stages were thus defined as the functions of daily temperature; such that an increase in temperature would result in an increase in the survival of all three adults stages and reduction in survival when temperature become lethal [22, 53]. The biological relationship between adult survival and temperature was assumed to be curvilinear [22, 53, 54].

Here, \({\varphi }_{0}, { \theta }_{0},\) and \({\alpha }_{0}\) refer to the baseline survival of unfed, blood-fed and gravid females respectively on a linear scale, under fixed temperature conditions of 27 ± 2 °C (insectary standard under which *An. funestus* have maximum survival [52, 54]), and assumes no blood meal limitation. The positive impact of temperature on all three life stages was represented by the coefficients \({\varphi }_{1},{ \theta }_{1}\) and \({\alpha }_{1}\) with an uninformative gamma prior with mean 12.5 and variance of 6.25. The coefficients \({ \varphi }_{2},{ \theta }_{2}\) and \({\alpha }_{2}\) correspond to the curvilinear effect of temperature on the survival of all three adult stages with their priors derived from the ratio between the linear coefficient and twice optimum temperature. This formulation ensured that the optimum temperature is fixed at (27 °C). The parameters \({\varepsilon }_{a}, { \varepsilon }_{f}, {\varepsilon }_{v}\) capture unexplained variation associated with survival during the distinct gonotrophic stages. These error terms (\({\varepsilon }^{*})\) were defined by normal priors with mean of 0 and a precision \({\sigma }^{*}\) from a gamma distribution with both shape and rate of 10 for unfed, blood-fed and gravid females.

#### Development between stages

The daily development probability from one life stage to the next was defined as the reciprocal of the development time (days) between the stages (assuming that all development times take longer than a day). An increase in temperature was assumed to reduce the development period of larvae [47, 53, 55, 56].

Specifically, \(L\left(t\right)\) is written as the function of temperature covariates:

Here \({C}_{0}\) corresponds to the baseline daily development period on a linear scale defined by an informative beta prior with range defined in Eq. 16–17 (Table 1). The coefficient \({C}_{1}\) explains the positive effect of temperature on larval development period, with its prior values derived from an uninformative gamma prior with mean 0.001 and standard deviation 0.001.

The development time for other life history stages (eggs and pupae) and the time between gonotrophic stages were assumed to be independent of temperature and other environmental covariates. The numbers of individuals (\({\mathcal{K}}_{m}\)) graduating from one stage to the next each day were modelled as a binomial process Eq. 11.

Here the rate \({r}_{m}\) is a development probability as defined in Eq. 11 for \(m\) stage, with assigned informative prior values as described through a generic prior in Eq. 16–17 (Table 1). Parameter \({\mathcal{W}}_{m-1}\left(t\right)\) refers to the number surviving the preceding life stage.

### Fecundity

The number of eggs laid at each time step was drawn from a Poisson distribution whose rate was the product of per-capita fecundity (number of eggs laid by blood-fed *An. funestus* under insectary conditions (\({b}_{0}\)), a penalized rate for the egg survival (\({s}_{e})\), the number of gravid mosquitoes (\({V}_{(t-1)})\) and ratio of females-males (assumed to be 0.5) as assessed at the pupae stage [19].

The error term \({\varepsilon }_{b}\) was defined by normal prior with mean of 0 and a precision from a gamma distribution with both shape and rate of 10.

### Observation-derived components of the Bayesian SSMs

Observations of the abundance of adults (unfed, bloodfed, gravid) \(A\) at time \(t\) were modelled as a normal distribution with varying daily means \(\mathrm{\overline{a} }\) determined by the biological model and a precision \(\tau\) representing observation error. A fixed coefficient of variation (\(\xi )\) for the daily observation process was assumed and assigned an uninformative prior with values between 0.1 and 0.9 (Table 1). The CDC light trap typically samples mosquitoes from populations of unknown size, for which the daily catch rates are difficult to quantify independently. A parameter \(\vartheta\) was therefore incorporated both into the precision \(\tau\) and daily varying means to account for an observed weekly periodicity in adult abundance, which was otherwise hard to interpret. This parameter was allowed to vary both by day of the week \(j\) and between the two populations \(k\) (2015 and 2018–19 datasets). The \(\vartheta\) values were derived from a logit function \(\mathrm{ exp}\left(\rho \right)/(1+\mathrm{exp}\left(\rho \right))\), with \(\rho\) defined from the uninformative normal prior with mean and standard deviation of 0 and 10 respectively. Therefore precision \(\tau\) can be written as \(\frac{1}{{({\xi }_{t}{\overline{a} }_{t}{\vartheta }_{jk})}^{2}}\) for all the adult stages. Thereafter, the observation abundance was estimated as follows:

#### Trap bias

The trapping method (CDC Light traps) primarily targets unfed host-seeking mosquitoes [57]. Blood fed and gravid mosquitoes are assumed to no longer host-seek, and represent a small proportion 0.5–3% of total collections of females caught [6, 22]). To account for these biases in sampling, a new parameter of “trap-biasness”\(\omega\) was added in the observation model for both precision \(\tau\) and varying daily means \({\overline{a} }_{\eta }\). The prior values for \(\omega\) were estimated from independent studies from the same locations [6, 58], and ranged from 0.05–0.15 (Table 1), with variations between the two life stages \(\eta\). Therefore, the observation model for blood-fed and gravid (\({A}_{\eta })\) was rewritten by modifying Eq. 14 as follows

#### Prior distributions

Since this model contains a large number of parameters, use of un-scaled informative priors restricted model convergence and mixing. A rescaled beta distributions of the informative priors [59] was opted and calibrated as follows:

where \(Y\) is a dummy variable that takes values in the interval [0, 1] with mean of 0.5 and standard deviation of 0.15, selected to provide low likelihood at the values 0 and 1. The values of \({X}_{min}\) and \({X}_{max}\) define the range of the parameter of interest as dictated by the prior information. Since the information on priors was provided in form of mean (\(\mu\)) and standard deviation (\(\sigma\)), the values were defined as \({X}_{min},{X}_{max}= \mu \pm 2\sigma\). Survival, development period, trap-biasness, variability in daily catches and fecundity parameters were all assigned priors according to Eq. 17.

### Model selection, model fitting and outputs

Model fitting was done using the R statistical software version 4.0.5 [60]. Population models were fitted using a Markov Chain and Monte Carlo sampling (MCMC) algorithm via the JAGS software [61] interfaced to R via the *runjags* package [62] (code provided in the Additional file 1: Fig. S1). To achieve convergence, the model with 6 chains was run in parallel for \({10}^{5}\) samples with a burn-in of \({10}^{5}\), keeping every 10th iteration for memory-saving reasons. Convergence was assessed by visual investigation of the trace plots, prior-posterior distribution using the *coda* package [63], effective sample sizes and the Gelman Rubin diagnostic [64]. Model comparisons were done using the deviance information criterion (DIC) [65], and the ones with the lowest DIC selected as the most preferred. The predicted and observed densities of *An. funestus* adult females were plotted to evaluate consistent prediction biases visually (Additional file 1: Fig. S2). Posterior means and 95% credible intervals for the key survival parameters, development period, density dependence, environmental covariates (temperature and rainfall) and fecundity were also reported to reveal different dynamical aspects of the system.

## Results

### Population trajectories and seasonal trends

Bayesian state-space model was used to describe the dynamics of wild populations of *An. funestus*. The full results, including summaries of posterior means for all the fitness and demographic parameters are reported in Table 1. The most parsimonious model (model-7, Table 3) included density dependence, and temperature and rainfall (current and lagged) impacts on larval survival, and the effect of temperature on larval development period. The only covariate that was not retained in the “model-7” was temperature impacts on adult survival. This model satisfactorily reconstructed the population dynamics of *An. funestu*s in the study villages, with all environmental covariate relaxation applied based on DIC selection. Population trajectories were estimated for all six *An. funestus* life history and gonotrophic stages after accounting for potential impacts of environmental covariates and density dependence (Fig. 3).

These trajectories reflect the annual trend in abundance spanning in periods from low or no rainfall to high rainfall. All trajectories show a relatively high abundance of *An. funestus* right after the rainy season, followed by reduced but sustained abundance during the dry period for all the life stages. After accounting for observation biases during sampling, the observed abundance of unfed, gravid and blood-fed groups largely falls within the credible intervals of the predicted values (Fig. 4).

### Survival and fecundity

Estimated *An. funestus* larval survival trajectories demonstrate substantial mean variability during the two seasons, with no clear pattern of seasonality (Fig. 5a and b), Table 1). Similarly, the survival trajectories of the adult stages (all gonotrophic states) were variable throughout the year, with daily survival rate ranging from 0.2 to 1.0 and not consistently differing between wet and dry seasons (Fig. 5c to h, Table 1). Per capita fecundity was estimated to be between 75 and 81 eggs per female *An. funestus* (Table 1). While the abundance of this species fluctuated seasonally, per capita fecundity remained consistent throughout the year (Fig. 5k and l).

Temperature was an important predictor of larval survival with a curvilinear relationship (ΔDIC = 138, Table 3, Fig. 6b), and that temperature has a positive monotonic relationship with larval development period (ΔDIC = 336, Table 3, Fig. 6a); with the larval development period estimated to last about 16 days on average. Additionally, daily rainfall was an important driver for the dynamics of *An. funestus* by reducing larval survival (ΔDIC = 5605, Table 3, Fig. 6c) with a negative monotonic relationship.

### Effects of density dependence

Density dependence was the only intrinsic feature incorporated in this dynamic model of *An. funestus*. The model was able to converge efficiently without crashing when density dependency was removed, suggesting this process plays a detectable but relatively minor role population regulation when compared with extrinsic factors (ΔDIC = 222, Table 3, Fig. 6e). To verify this, a simulation was run and discovered that the estimated density dependence was actually quite low, even when simulating with a single dataset. The model fitting process also suggested the interaction parameter (\({\beta }_{3})\) between larvae density \((D)\) and one week cumulative rainfall \((Q)\) contributes to *An. funestus* dynamics by positively increasing larval survival (ΔDIC = 38, Table 3, Fig. 6d).

## Discussion

Advanced methods used in mainstream ecological studies was adapted and fitted a state-space model (SSM) to field and laboratory data to accurately reconstruct population dynamics of wild population of *An. funestus.* The SSM inferred the trajectories of multiple life-cycle and gonotrophic stages of wild *An. funestus* females. This allows the reconstruction of the observed trajectories of larvae and adult females for the wild *An. funestus* in Tanzania for the first time. This analysis indicated that the dynamics of *An. funestus* were best explained in a model that included density dependency, temperature (curvilinear relationship) and rainfall (negative monotonic relationship) on larval survival, and temperature on the larval development period (positive monotonic relationship). In contrast, model fit was not improved by incorporating temperature dependency into adult survival (all gonotrophic stages). *Anopheles funestus* abundance vary seasonally between wet and dry but demographics rates (i.e. survivals, fecundity and development period) did not vary after accounting for the impact of environmental covariates and density dependence. These results are very useful for generating hypotheses about the nature and relative magnitude of drivers of *An. funestus* population dynamics in the wild. This model can be extended to include a component on malaria dynamics in humans; or to compare the efficacy and effectiveness of different interventions in combination or singly. This would allow more sophisticated evaluation of suitability of *An. funestus*-specific interventions; including prediction of the potential combined effect of strategies that acting at different life-cycle stages and/or target different demographic processes (e.g. survival versus fecundity).

Extrinsic covariates such as rainfall and temperature were all hypothesised to be the main drivers for the dynamics of this vector species. This study supports the hypothesis that rainfall is a significant driver of the population dynamics of wild *An. funestus*. Overall, the abundance of all life stages were relatively higher in rainy compared to dry periods of the year as previously documented [22, 66,67,68]. Rainfall covariates were directly included in the larval survival model since it is the only stage on which rainfall was hypothesized to have a significant impact. Daily larval survival as estimated by the SSM showed high variability both within seasons and across the year. There was support for a monotonic association between rainfall and larval survival; characterized as reduction in larval survival during periods of heavy rainfall [69]. Time lags have been used to assess rainfall impacts on the dynamics of these vectors [22]. *Anopheles funestus* abundance have been shown to be positively associated with the cumulative lag rainfall [22]. Here one week cumulative rainfall was included in the model to account for its effect on survival. Similar to other vectors of malaria transmission such as *An. gambiae*, rainfall have always been considered as the main factor regulating the dynamics, despite ecological differences between the two vector species [22, 46, 70, 71].

The SSM also provided support for the hypothesis that temperature is an important driver of *An. funestus* dynamics; although the nature of temperature effects was complex and variable between life history stages. For example, temperature was associated with both larval survival and development, but not adult survival or fitness. Furthermore the estimated impacts of temperature on larval ecology were complex; with the SSM suggesting a curvilinear relationship with survival but a positive monotonic impact with the larval development period. These findings validate the prior studies that demonstrated that temperature had a curvilinear influence on *Anopheles* larval survival, with a rise in temperature above/below the optimum lowering survival [47, 48, 72]. The larval period of *An. gambiae* is temperature dependent [55, 72]; thus this model incorporated a positive monotonic effect such that development is fastest when temperature is high and just below maximum threshold for larval development [72, 73]. In the final model the effect of temperature on gonotrophic stages was not found to be an important driver for the dynamic of this species thus left out during model selection process.

Little is known about the effect of density dependence on *An. funestus* due to its ecology and reliance on the large semi-permanent and permanent breeding habitats [20, 21]. However, density dependence is already well-known to be an important driver for dynamic of other malaria vectors like *An. gambiae* [25, 27,28,29, 31, 67, 74] and other non-malaria vectors like *Aedes aegypti* [75, 76]. Variations in densities during the aquatic stages of *An. gambiae s.l.* have been found to affect adult fitness [28, 77]. The SSM fit better when density dependence of larval survival was included, though the relative magnitude of this process was quite small and likely to have minor impact on the overall dynamics of *An. funestus* populations (Fig. 6e). These findings suggest that *An. funestus* populations are likely to be regulated more by extrinsic than intrinsic processes. These findings corroborate the original hypotheses about density dependence having a weaker regulatory role in this species on the basis of the types of larval habitats (i.e. larger and more permanent habitat [20, 78]), which can likely sustain higher resources and thus reduce competition than in *An. gambiae*. This is the first report documenting the role of density dependence on the dynamics of the wild populations of *An. funestus*. Now that colonies are becoming more feasible, more thorough investigation on the role of density dependence in the dynamic of *An. funestus* is prerequisite.

In addition to highlighting potential drivers of *An. funestus* populations, the SSM here generated plausible estimates of key demographic and life-history process in the wild. This model estimated that *An. funestus* larvae takes an average of 15.6–16.1 days to grow from first instar larvae to pupae; which is relatively long compared to the other major vectors in the *An. gambiae* complex (9–11 days [48, 55]). This apparently longer development period of *An. funestus* may be a product of their adaptation to more permanent, year-round breeding habitats that are unlikely to dry up; thus reducing selection for rapid development. The SSM estimated that the daily survival of wild *An. funestus* larvae could be as high as 0.95, compared to the 0.83 [0.80, 0.86] mean daily survival rate of the known vector of malaria transmission *An. gambiae* [48, 79]. This matches observations from insectary experiments in which *An. funestus* larvae have higher survival than *An. gambiae* [19, 38, 80]. Given the apparently higher rates of survival in *An. funestus* than in *An. gambiae*, these findings suggest that more lethal intervention may be required to control *An. funestus* both at larvae and adults stages.

The impact of any vector control largely depends on the ecology of the specific vector species. Differences in ecology between *An. gambiae* and *An. funestus* are likely to affect the relative impact of interventions. For instance, *An. gambiae* prefer breeding in small and temporary habitats which dry up quickly when there is no rainfall which is opposite to *An. funestus* habitats. Despite the fact that *An. funestus* habitats are "few, fixed, and findable" and might be easily targeted for larviciding [20] during the dry season, treating habitats such as rivers or bigger ponds could pose logistical challenges. The persistence of *An. funestus* throughout the year even during the driest periods suggest this vector is less seasonal compared to *An. gambiae s.l*., which experience much more dramatic “boom and bust” dynamics in relation to seasonal rains [6, 22, 68]. If interpreted together with the observation that survival estimates were not seasonal, the model suggests that this species is likely responsible for year-round malaria transmission throughout the year, while other species, which mostly occupy temporary habitats may be responsible for any apparent seasonality in transmission.

Models of vector population dynamics can provide a useful guide for the selection of optimal vector control strategies; particular through enabling more focal investigation of the benefits of seasonal or spatial targeting and use of combined versus single interventions. Despite its complexity, this population dynamics model provides a useful framework for investigation of the stability of *An. funestus* populations. With additional data, this model can be further refined to include additional modifications related to vector ecology and behaviour that may impact intervention (e.g. host choice and its impacts on fitness, predation during larval or adult phase and spatial components). Such further elucidation may increase the predictive accuracy of this SSM in specific contexts, but even the more general framework developed here have flexibility to introduce stage-specific mortality effects expected from different types of vector control [16, 32, 33]. For example, this framework could be used to model the impact of combined interventions including those that target adult females (insecticide-treated nets (ITNs), IRS) and larviciding; and assessment of how mortality varies with different coverage [16, 32]. It can also be used to investigate the possible response of vector population to climate change anticipated in Tanzania and other African countries. An important limitation of this study is lack of knowledge on what percentage of *An. funestus* mosquito population is sampled by the trap, which is important for understanding the relative magnitude of demographic stochasticity in modelled dynamics. This highlights the need to explicitly incorporate this source of uncertainty into vector and transmission dynamics; including the need for further calibration and standardization of the efficiency and biases associated with particular mosquito trapping methods.

## Conclusions

This study used Bayesian State Space Models (SSM) parameterized with empirical data to quantify key demographic and fitness processes underpinning the population dynamics of *An. funestus* in Tanzania. This is the first use of SSM to understand the population dynamic of the wild vector of residual malaria transmission, *An. funestus* in Tanzania. The model structure allowed investigation of the relative importance of seasonally-varying environmental covariates (i.e. rainfall and temperature) and density dependence; providing some support for both processes although the magnitude of the former was much greater than the latter. The ability of this model to accurately reconstruct the seasonal dynamics and demography of *An. funestus* indicate its value for simulating the response of these populations to vectors control measures applied either individually or in combination. Additionally, the relatively limited evidence of seasonality in key fitness and demographic rates further corroborate evidence that this vector species can facilitate efficient year-round transmission of malaria. Finally, this model also highlights the clear importance of accounting for regional and daily observation biases when modelling mosquito population dynamics.

## Availability of data and materials

The R-codes used for model development are freely available in the supplementary file. The time series data is available upon reasonable request to the author.

## References

Gillies M, Meillon D (1968) The Anophelinae of Africa south of the Sahara (Ethiopian zoogeographical region). Johannesburg: South Afr Inst Med Res. 1968. p. 343.

Sinka ME, Bangs MJ, Manguin S, Rubio-Palis Y, Chareonviriyaphap T, Coetzee M, et al. A global map of dominant malaria vectors. Parasit Vectors. 2012;5:69.

Russell TL, Govella NJ, Azizi S, Drakeley CJ, Kachur SP, Killeen GF, et al. Increased proportions of outdoor feeding among residual malaria vector populations following increased use of insecticide-treated nets in rural Tanzania. Malar J. 2011;10:80.

Hunt RH, Brooke BD, Pillay C, Koekemoer LL, Coetzee M. Laboratory selection for and characteristics of pyrethroid resistance in the malaria vector

*Anopheles funestus*. Med Vet Entomol. 2005;19:271–5.Braack LE, Coetzee M, Hunt RH, Biggs H, Cornel A, Gericke A. Biting pattern and host-seeking behavior of

*Anopheles arabiensis*(Diptera: Culicidae) in northeastern South Africa. J Med Entomol. 1994;31:333–9.Kaindoa EW, Matowo NS, Ngowo HS, Mkandawile G, Mmbando A, Finda M, et al. Interventions that effectively target

*Anopheles funestus*mosquitoes could significantly improve control of persistent malaria transmission in south-eastern Tanzania. PLoS ONE. 2017;12: e0177807.Lwetoijera D, Harris C, Kiware SS, Dongus S, Devine GJ, McCall PJ, et al. Increasing role of

*Anopheles funestus*and*Anopheles arabiensis*in malaria transmission in the Kilombero Valley, Tanzania. Malar J. 2014;13:331.Garrett-Jones C, Shidrawi GR. Malaria vectorial capacity of a population of

*Anopheles gambiae*: an exercise in epidemiological entomology. Bull World Health Organ. 1969;40:531–45.Seyoum A, Sikaala CH, Chanda J, Chinula D, Ntamatungiro AJ, Hawela M, et al. Human exposure to anopheline mosquitoes occurs primarily indoors, even for users of insecticide-treated nets in Luangwa Valley. Southeast Zambia Parasit Vectors. 2012;5:101.

Riveron JM, Chiumia M, Menze BD, Barnes KG, Irving H, Ibrahim SS, et al. Rise of multiple insecticide resistance in

*Anopheles funestus*in Malawi: a major concern for malaria vector control. Malar J. 2015;14:344.Hunt RH, Edwardes M, Coetzee M. Pyrethroid resistance in southern African An

*opheles funestus*extends to Likoma Island in Lake Malawi. Parasit Vectors. 2010;3:122.McCann RS, Ochomo E, Bayoh MN, Vulule JM, Hamel MJ, Gimnig JE, et al. Reemergence of

*Anopheles funestus*as a vector of*Plasmodium falciparum*in Western Kenya after long-term implementation of insecticide-treated bed nets. Am J Trop Med Hyg. 2014;90:597–604.Bayoh MN, Mathias DK, Odiere MR, Mutuku FM, Kamau L, Gimnig JE, et al.

*Anopheles gambiae*: historical population decline associated with regional distribution of insecticide-treated bed nets in western Nyanza Province, Kenya. Malar J. 2010;9:62.Killeen GF, Seyoum A, Sikaala C, Zomboko AS, Gimnig JE, Govella NJ, et al. Eliminating malaria vectors. Parasit Vectors. 2013;6:172.

Okumu F, Finda M. Key characteristics of residual malaria transmission in two districts in South-Eastern Tanzania-implications for improved control. J Infect Dis. 2021;223:S143–54.

Fillinger U, Kannady K, William G, Vanek MJ, Dongus S, Nyika D, et al. A tool box for operational mosquito larval control: preliminary results and early lessons from the Urban Malaria Control Programme in Dar es Salaam. Tanzania Malar J. 2008;7:20.

Killeen GF, Chitnis N, Moore SJ, Okumu FO. Target product profile choices for intra-domiciliary malaria vector control pesticide products: repel or kill? Malar J. 2011;10:207.

Griffin JT, Hollingsworth TD, Okell LC, Churcher TS, White M, Hinsley W, et al. Reducing

*Plasmodium falciparum*malaria transmission in Africa: a model-based evaluation of intervention strategies. PLoS Med. 2010;7: e1000324.Ngowo HS, Hape EE, Matthiopoulos J, Okumu FO. Fitness characteristics of the malaria vector,

*Anopheles funestus*, during an attempted laboratory colonization. Malar J. 2020;20:148.Nambunga IH, Ngowo HS, Mapua SA, Hape EE, Msugupakulya BJ, Msaky DS, et al. Aquatic habitats of the malaria vector

*Anopheles funestus*in rural south-eastern Tanzania. Malar J. 2020;19:219.Gimnig JE, Ombok M, Kamau L, Hawley WA. Characteristics of larval Anopheline (Diptera: Culicidae) habitats in Western Kenya. J Med Entomol. 2001;38:282–8.

Ngowo HS, Kaindoa EW, Matthiopoulos J, Ferguson HM, Okumu FO. Variations in household microclimate affect outdoor-biting behaviour of malaria vectors. Wellcome Open Res. 2017;2:102.

Minakawa N, Sonye G, Mogi M, Yan G. Habitat characteristics of

*Anopheles gambiae s.s*. larvae in a Kenyan highland. Med Vet Entomol. 2004;18:301–5.Kabbale FG, Akol AM, Kaddu JB, Onapa AW. Biting patterns and seasonality of

*Anopheles gambiae sensu lato*and*Anopheles funestus*mosquitoes in Kamuli District, Uganda. Parasit Vectors. 2013;6:340.Russell TL, Lwetoijera DW, Knols BGJ, Takken W, Killeen GF, Ferguson HM. Linking individual phenotype to density-dependent population growth: the influence of body size on the population dynamics of malaria vectors. Proc R Soc B Biol Sci. 2011;278:3142–51.

Nowicki P, Bonelli S, Barbero F, Balletto E. Relative importance of density-dependent regulation and environmental stochasticity for butterfly population dynamics. Oecologia. 2009;161:227–39.

Gimnig JE, Ombok M, Otieno S, Michael G, Vulule JM, Walker ED, et al. Density-dependent development of

*Anopheles gambiae*(Diptera : Culicidae) larvae in artificial habitats. J Med Entomol. 2002;39:162–72.Ng’habi KB, John B, Nkwengulila G, Knols BGJ, Killeen GF, Ferguson HM. Effect of larval crowding on mating competitiveness of

*Anopheles gambiae*mosquitoes. Malar J. 2005;4:49.Koenraadt CJM, Majambere S, Hemerik L, Takken W. The effects of food and space on the occurrence of cannibalism and predation among larvae of

*Anopheles gambiae s.l*. Entomol Exp Appl. 2004;112:125–34.Charlwood JD. May the force be with you: measuring mosquito fitness in the field. Ecol Asp Appl Genet Modif Mosq. 2003;2:47–62.

Henderson PA, Magurran AE. Direct evidence that density-dependent regulation underpins the temporal stability of abundant species in a diverse animal community. Proc R Soc B Biol Sci. 2014;281:20141336.

Viana M, Ng’habi K, Lyimo I, Ferguson HM, Matthiopoulos J, Killeen G. Mesocosm experiments reveal the impact of mosquito control measures on malaria vector life history and population dynamics. Sci Rep. 2018;8:13949.

Viana M, Hughes A, Matthiopoulos J, Ranson H, Ferguson HM. Delayed mortality effects cut the malaria transmission potential of insecticide-resistant mosquitoes. Proc Natl Acad Sci USA. 2016;113:8975–80.

Lyimo IN, Haydon DT, Russell TL, Mbina KF, Daraja AA, Mbehela EM, et al. The impact of host species and vector control measures on the fitness of African malaria vectors. Proc Biol Sci. 2013;280:20122823.

Mweresa CK, Omusula P, Otieno B, van Loon JJA, Takken W, Mukabana WR. Molasses as a source of carbon dioxide for attracting the malaria mosquitoes

*Anopheles gambiae*and*Anopheles funestus*. Malar J. 2014;13:160.Lyimo IN, Haydon DT, Mbina KF, Daraja AA, Mbehela EM, Reeve R, et al. The fitness of African malaria vectors in the presence and limitation of host behaviour. Malar J. 2012;11:425.

Ferguson HM, Ng’habi KR, Walder T, Kadungula D, Moore SJ, Lyimo I, et al. Establishment of a large semi-field system for experimental study of African malaria vector ecology and control in Tanzania. Malar J. 2008;7:158.

Okoye PN, Brooke BD, Hunt RH, Coetzee M. Relative developmental and reproductive fitness associated with pyrethroid resistance in the major southern African malaria vector, Anopheles funestus. Bull Entomol Res. 2007;97:599–605.

Newman KB, Buckland ST, Lindley ST, Thomas L, Fernández C. Hidden process models for animal population dynamics. Ecol Appl. 2006;16:74–86.

Millar R. Bayesian state-space modeling of age-structured data: fitting a model is just the beginning. Can J Fish Aquat Sci. 2000;50:43–50.

Viana M, Cleaveland S, Matthiopoulos J, Halliday J, Packer C, Craft ME, et al. Dynamics of a morbillivirus at the domestic–wildlife interface: canine distemper virus in domestic dogs and lions. Proc Natl Acad Sci USA. 2015;112:1464–9.

Rivot E, Prévost E, Parent E, Baglinière JL. A Bayesian state-space modelling framework for fitting a salmon stage-structured population dynamic model to multiple time series of field data. Ecol Modell. 2004;179:463–85.

Mboera LEG. Sampling techniques for adult Afrotropical malaria vectors and their reliability in the estimation of entomological inoculation rate. Tanzan Health Res Bull. 2005;7:117–24.

Gillies MT. The duration of the gonotrophic cycle in

*Anopheles gambiae*and*Anopheles funestus*, with a note on the efficiency of hand catching. East Afr Med J. 1953;30:129–35.Rúa GL, Quiñones ML, Vélez ID, Zuluaga JS, Rojas W, Poveda G, et al. Laboratory estimation of the effects of increasing temperatures on the duration of gonotrophic cycle of

*Anopheles albimanus*(Diptera: Culicidae). Mem Inst Oswaldo Cruz. 2005;100:515–20.Fillinger U, Sonye G, Killeen GF, Knols BGJ, Becker N. The practical importance of permanent and semipermanent habitats for controlling aquatic stages of

*Anopheles gambiae*sensu lato mosquitoes: Operational observations from a rural town in western Kenya. Trop Med Int Health. 2004;9:1274–89.Bayoh MN, Lindsay SW. Temperature-related duration of aquatic stages of the Afrotropical malaria vector mosquito

*Anopheles gambiae*in the laboratory. Med Vet Entomol. 2004;18:174–9.Lyimo EO, Takken W, Koella JC. Effect of rearing temperature and larval density on larval survival, age at pupation and adult size of

*Anopheles gambiae*. Entomol Exp Appl. 1992;63:265–71.Beier JC, Copeland R, Oyaro C, Masinya A, Odago WO, Oduor S, et al.

*Anopheles gambiae*complex egg-stage survival in dry soil from larval development sites in western Kenya. J Am Mosq Control Assoc. 1990;6:105–9.Koenraadt CJM, Majambere S, Hemerik L, Takken W. Cannibalism and predation among larvae of

*Anopheles gambiae s.l*. Entomol Exp Appl. 2004;112:125–34.Paaijmans KP, Heinig RL, Seliga RA, Blanford JI, Blanford S, Murdock CC, et al. Temperature variation makes ectotherms more sensitive to climate change. Glob Chang Biol. 2013;19:2373–80.

Paaijmans KP, Imbahale SS, Thomas MB, Takken W. Relevant microclimate for determining the development rate of malaria mosquitoes and possible implications of climate change. Malar J. 2010;9:196.

Blanford JI, Blanford S, Crane RG, Mann ME, Paaijmans KP, Schreiber KV, et al. Implications of temperature variation for malaria parasite development across Africa. Sci Rep. 2013;3:1300.

Paaijmans KP, Thomas MB. The influence of mosquito resting behaviour and associated microclimate for malaria risk. Malar J. 2011;10:183.

Kirby MJ, Lindsay SW. Effect of temperature and inter-specific competition on the development and survival of

*Anopheles gambiae sensu stricto*and*An. arabiensis*larvae. Acta Trop. 2009;109:118–23.Paaijmans KP, Blanford S, Bell AS, Blanford JI, Read AF, Thomas MB. Influence of climate on malaria transmission depends on daily temperature variation. Proc Natl Acad Sci USA. 2010;107:15135–9.

Mboera L. Sampling techniques for adult Afrotropical malaria vectors and their reliability in the estimation of entomological inoculation rate. Tanzan J Health Res. 2006;7:117–24.

Msugupakulya BJ, Kaindoa EW, Ngowo HS, Kihonda JM, Kahamba NF, Msaky DS, et al. Preferred resting surfaces of dominant malaria vectors inside different house types in rural south-eastern Tanzania. Malar J. 2020;19:22.

Gelman A, Jakulin A, Pittau MG, Su YS. A weakly informative default prior distribution for logistic and other regression models. Ann Appl Stat. 2008;2:1360–83.

R Development Core Team. R: A language and environment for statistical computing. R Found Stat Comput. 2021.

Plummer M. JAGS : a program for analysis of Bayesian graphical models using Gibbs Sampling JAGS : Just Another Gibbs Sampler. 2003.

Denwood MJ. runjags : an R package providing interface utilities, model templates, parallel computing methods and additional distributions for MCMC models in JAGS. J Stat Softw. 2016;71:1–25.

Plummer M, Best N, Cowles K, Vines K. CODA: Convergence Diagnosis and Output Analysis for MCMC. R News. 2006.

Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences. Stat Sci. 1992. https://doi.org/10.1214/ss/1177011136.

Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A. Bayesian measures of model complexity and fit. J R Stat Soc Ser C Appl Stat. 2002;64:583–639.

Kreppel KS, Viana M, Main BJ, Johnson PCD, Govella NJ, Lee Y, et al. Emergence of behavioural avoidance strategies of malaria vectors in areas of high LLIN coverage in Tanzania. Sci Rep. 2020;10:14527.

Smith T, Charlwood JD, Takken W, Tanner M, Spiegelhalter DJ. Mapping the densities of malaria vectors within a single village. Acta Trop. 1995;59:1–18.

Charlwood JD, Vij R, Billingsley PF. Dry season refugia of malaria-transmitting mosquitoes in a dry savannah zone of east Africa. Am J Trop Med Hyg. 2000;62:726–32.

Kelly-Hope AL, Hemingway J, McKenzie FE. Environmental factors associated with the malaria vectors

*Anopheles gambiae*and*Anopheles funestus*in Kenya. Malar J. 2009;8:268.Lyons CL, Coetzee M, Terblanche JS, Chown SL. Thermal limits of wild and laboratory strains of two African malaria vector species,

*Anopheles arabiensis*and*Anopheles funestus*. Malar J. 2012;11:226.Zhou G, Munga S, Minakawa N, Githeko AK, Yan G. Spatial relationship between adult malaria vector abundance and environmental factors in western Kenya highlands. Am J Trop Med Hyg. 2007;77:29–35.

Paaijmans KP, Huijben S, Githeko AK, Takken W. Competitive interactions between larvae of the malaria mosquitoes

*Anopheles arabiensis*and*Anopheles gambiae*under semi-field conditions in western Kenya. Acta Trop. 2009;109:124–30.Lyons CL, Coetzee M, Chown SL. Stable and fluctuating temperature effects on the development rate and survival of two malaria vectors,

*Anopheles arabiensis*and*Anopheles funestus*. Parasit Vectors. 2013;6:104.Muriu SM, Coulson T, Mbogo CM, Godfray HCJ. Larval density dependence in

*Anopheles gambiae s.s*., the major African vector of malaria. J Anim Ecol. 2013;82:166–74.Yang GJ, Brook BW, Whelan PI, Cleland S, Bradshaw CJA. Endogenous and exogenous factors controlling temporal abundance patterns of tropical mosquitoes. Ecol Appl. 2008;18:2028–40.

Hancock PA, White VL, Callahan AG, Godfray CHJ, Hoffmann AA, Ritchie SA. Density-dependent population dynamics in

*Aedes**aegypti*slow the spread of wMel Wolbachia. J Appl Ecol. 2016;53:785–93.Porretta D, Mastrantonio V, Crasta G, Bellini R, Comandatore F, Rossi P, et al. Intra-instar larval cannibalism in

*Anopheles**gambiae*(s.s.) and*Anopheles**stephensi*(Diptera: Culicidae). Parasit Vectors. 2016;9:566.Tuno N, Githeko A, Yan G, Takagi M. Interspecific variation in diving activity among

*Anopheles gambiae*Giles,*An. arabiensis*Patton, and*An. funestus*Giles (Diptera: Culicidae) larvae. J Vector Ecol. 2007;32:112–7.Matthews J, Bethel A, Osei G. An overview of malarial

*Anopheles*mosquito survival estimates in relation to methodology. Parasit Vectors. 2020;13:233.Takken W, Charlwood JD, Billingsley PF, Gort G. Dispersal and survival of

*Anopheles**funestus*and*A. gambiae s.l.*(Diptera: Culicidae) during the rainy season in southeast Tanzania. Bull Entomol Res. 1998;88:561.

## Acknowledgements

We are grateful to member of IHI-VectorSphere insectary’s Joseph Mgando, Dickson, Mwasheshi, Godfrey Matanila, Neema Nombo, Rukiyah Mohammed, Johnson Kyeba and Khamis Bwanally for their contribution in different occasion during mosquito rearing, colonization and collection of wild population data. The earlier version of the model development benefited from the contributions of Luca Nelli, Fergus Chadwick, Heather McDevitt, Crinan Jarrett, Yacob Haddou and Prashanth Selvaraj. We also thank the three anonymous reviewers for their contributions.

## Funding

This work was supported by Howard Hughes Medical Institute (HHMI)-Gates Foundation (Grants: OPP1099295) and Gates Foundation (Grants: INV-002138) both awarded to FOO and HMF which support HSN.

## Author information

### Authors and Affiliations

### Contributions

HSN, FOO, HMF and JM designed the study. HSN, IHM and EEH conducted the lab experiment and collected data used in model development. HSN and JM conducted the formal model development, fitting and validation. HSN drafted the manuscript. FOO, HMF and JM reviewed the manuscript and provide supervision. All authors read and approved the final manuscript.

### Corresponding author

## Ethics declarations

### Ethics approval and consent to participate

This study was approved by Ifakara Health Institute Review Board (Ref: IHI/IRB/No: 007-2018) and the Medical Research Coordinating Committee (MRCC) at the National Institute for Medical Research-NIMR (Ref: NIMR/HQ/R.8a/Vol.IX/2895). Individual verbal and written consent was also obtained from household owner where CDC light traps were places for collecting mosquitoes.

### Consent for publication

Permission to publish was obtained from Medical Research Coordinating Committee (MRCC) at the National Institute for Medical Research (Ref. No: NIMR/HQ/P.12 VOL XXXIII/116).

### Competing interests

The authors declare no conflicts of interest.

## Additional information

### Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supplementary Information

**Additional file 1: Figure S1.**

SSMs model development R codes. **Figure S2.** Goodness-of-fit: Observed versus predicted unfed, bloodfed and gravid densities across all populations. Adjusted R-squared, intercept and slope values are from a linear model of the predicted against observed values. Dotted lines correspond to 1:1 line. Left column (a,c,d) is data collected from June 2018 to May 2019 and right column (b,d,e) is data from Jan-Dec 2015. Grey area is the period with rainfall. **Figure S3.1.** Prior (orange histogram) and posterior (blue histogram) distribution of the main baseline and observational parameters in the state-space model. **Figure S3.2.** Prior (orange histogram) and posterior (blue histogram) distribution of the main environmental covariates parameters in the state-space model.

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

## About this article

### Cite this article

Ngowo, H.S., Okumu, F.O., Hape, E.E. *et al.* Using Bayesian state-space models to understand the population dynamics of the dominant malaria vector, *Anopheles funestus* in rural Tanzania.
*Malar J* **21**, 161 (2022). https://doi.org/10.1186/s12936-022-04189-4

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s12936-022-04189-4

### Keywords

- Anopheles
*funestus* - State space model
- Population dynamic
- Seasonality
- Abundance
- Density dependence