Model-based analysis of experimental data from interconnected, row-configured huts elucidates multifaceted effects of a volatile chemical on Aedes aegypti mosquitoes

Abstract

Background

Insecticides used against Aedes aegypti and other disease vectors can elicit a multitude of dose-dependent effects on behavioral and bionomic traits. Estimating the potential epidemiological impact of a product requires thorough understanding of these effects and their interplay at different dosages. Volatile spatial repellent (SR) products come with an additional layer of complexity due to the potential for altered movement of affected mosquitoes and diffusion of volatile particles of the product beyond the treated house. Here, we propose a paired experimental design and statistical inference framework for estimating these nuanced effects of volatile SRs.

Method

We fitted a continuous-time Markov chain model in a Bayesian framework to data on marked mosquitoes released in interconnected experimental huts conducted in Iquitos, Peru. We estimated the effects of two dosages of transfluthrin on Ae. aegypti behaviors associated with human-vector contact: repellency, exiting and knockdown in the treated space and in connected, adjacent huts. We validated the framework using simulated data.

Results

The odds of a female Ae. aegypti being repelled, and thus prevented from entering a treated hut (H_T), increased at both dosages (low dosage: odds = 1.64, 95% highest density interval (HDI) = 1.30–2.09; high dosage: odds = 1.35, HDI = 1.04–1.67). The relative risk of exiting from the treated hut was reduced (low: RR = 0.70, HDI = 0.62–1.09; high: RR = 0.70, HDI = 0.40–1.06), with this effect carrying over to untreated spaces two huts away from the treated hut (H₂) (low: RR = 0.79, HDI = 0.59–1.01; high: RR = 0.66, HDI = 0.50–0.87). Knockdown rates were increased in both treated and downstream huts, particularly under high dosage (H_T: RR = 8.37, HDI = 2.11–17.35; H₁: RR = 1.39, HDI = 0.52–2.69; H₂: RR = 2.22, HDI = 0.96–3.86).

Conclusions

Our framework is effective at elucidating multiple effects of volatile chemicals used in SR products, as well as their downstream effects. For the examined formulations of transfluthrin, we found notable dose-dependent effects on repellency, movement and knockdown that carry over to adjacent, untreated spaces.

Background

Insecticidal strategies against adult mosquitoes have been used extensively in the control of mosquito-borne diseases [1]. However, certain mosquito behaviors, such as outdoor and daytime biting, challenge the efficacy of traditional control tools like insecticide treated nets (ITNs) and indoor residual spraying (IRS) [2]. The evolution of physiological resistance to insecticides [3] and behavioral adaptation of mosquitoes [4, 5] also pose limitations to the effectiveness of such products.

The effect of vector control products often goes beyond their acute lethal effects. For example, ITNs can elicit knockdown with potential for mosquito recovery and can divert mosquitoes away from a protected human to alternate hosts [6,7,8]. Volatile chemicals such as transfluthrin and metofluthrin can be delivered in high dosages and result in high lethality or can be formulated at lower dosages where acute toxicity is attenuated and other sublethal effects, such as incapacitated host-attraction and/or blood-feeding, are elicited instead, as was described previously for residual pyrethroids [9, 10]. Currently, the term “spatial repellency” is used to describe a range of behaviors that products with volatile chemicals, including spatial repellents (SR), may invoke [11], including repellency (reduced entry), diversion (product-induced movement to an untreated space), knockdown, mortality, irritancy (increase in exiting), attraction-inhibition and reduced biting [12,13,14]. These modes of action can have a concerted impact on disease transmission on an individual and community level [15,16,17,18].

Mark-release-recapture-type (MRR) experimental hut studies, in which marked mosquitoes are released in or outside an experimental hut system and monitored using traps and knockdown rates, offer unique opportunities to elucidate dosage and behavioral effects of SRs by measuring lethality, repellency and irritancy of a target vector species [10, 19,20,21]. However, studies such as these have not yet provided the granularity required to disentangle distinct behavioral and bionomic effects. The primary challenge associated with the design and interpretation of these studies is that each mosquito is only observed once: when knocked down or when trapped in entry or exit traps. This leaves movement trajectories in between release and recapture locations unobserved, making it challenging to quantify the relative contributions of multiple competing effects that could account for observed individual-level outcomes under a multitude of equally plausible scenarios. One recent study [22] showed that even short periods of transient exposure to volatile SRs can have significant, and sometimes delayed, effects on vectors. Such unobserved effects may compromise traditional statistical analyses.

Models used for the analysis of MRR data have a long history in ecology [23,24,25,26,27]. Originally developed to estimate survival probabilities and population sizes [28], they are now increasingly being used to inform spatial processes [29]. These models partition animal movement trajectories into states (e.g. breeding or foraging), with multi-state MRR models accounting for the probability of the animal occupying any of the possible states at a given time. Given sufficient information from sampling at multiple points in time and appropriate model constraints, these models can be extended for parameter estimation in the presence of unobserved states [30]. Bayesian methods are increasingly being applied to these types of problems given their treatment of all quantities as random variables [31,32,33]. These methods allow for formal treatment and quantification of parameter uncertainty and they allow researchers to explicitly build on previous studies.

Here, we present a new approach to infer effects of SR products on adult female Ae. aegypti by developing a hierarchical Bayesian model and applying it to a MRR study uniquely designed to measure SR product-induced effects resulting from a single treated home amongst a series of adjacent, untreated homes. The design mimicking an anticipated SR operational scenario in a typical semi-urban, dengue-endemic environment whereby product is implemented at < 100% household coverage and houses are adjoined, thereby forming shared interior spaces. This analytical approach was necessitated by the fact that the MRR study was conducted in a set of five interconnected huts, meaning that the time spent by individual mosquitoes in each of the five huts could not be directly observed. We first demonstrate the accuracy of this approach using data simulated under the same design as in our field experiments. We then demonstrate the dose-dependency of knockdown, repellency and exiting effects of the SR in both treated and untreated huts. We discuss the potential use of this framework to inform the projected impact and implementation of SRs and other vector control tools with volatile chemicals.

Methods

The study was performed at the Instituto Veterinario de Investigaciones in Iquitos, Peru (73.2°W, 7.3°S) during May-June, 2011.

Product

Technical grade transfluthrin (Sigma-Aldrich 105 Co. LLC, St. Louis, MO, USA), a volatile pyrethroid insecticide, was applied to cotton strips at 1/16th (0.0025 g/m²) and 1/8th (0.005 g/m²) dilutions of a field application rate (0.04 g/m²) using previously established protocols [34]. Control strips of matched cotton material were treated with acetone alone. Cotton material was applied to the interior walls of a designated experimental hut at 25% surface area coverage using magnets and metal frames [34].

Experimental huts

A unique experimental hut configuration was used for this study. Five independent structures were positioned adjacent to one another in a single row creating adjoining walls (Figs. 1, 2). Side eave gaps of each structure were open, subsequently allowing a continuum of indoor space among all five huts available for mosquito and volatile chemical movement. This design was deliberate to mimic the housing configuration common in the study location in Iquitos, Peru and other dengue-endemic areas (i.e. semi-urban settings in resource-poor, tropical areas). The huts measured 4 × 6 m and had 2 m high sidewalls. Each sidewall had an eave space 40 cm high and 6 m long (one-fifth of wall surface) (Fig. 2). The eave passage between adjoining huts was fitted with netting on all sides except those facing the interior hut space. In essence, mosquitoes could move between hut interiors but not enter the exterior space between huts. Each hut had two windows (one each on the front and back walls) equipped with exit traps and each hut had two doors (one of each on the front and back walls). All window and door openings were equipped with exit intervention traps. The two outermost huts had additional eave exit intervention traps on sidewalls (Fig. 2). Hut construction materials and structural design were based on previous MRR hut studies [19, 20].

Fig. 1

Experimental huts in row house configuration. Iquitos, Peru

Fig. 2

Longitudinal (a) and transversal (b) images of the interconnected huts with door, window and eave traps (c-e) and the eaves between adjoining walls (f)

Mosquitoes

Female Ae. aegypti test populations (F₁-F₂ generations) of 5–7 days-old were reared from field-collected larvae following previously established protocols [35]. Mosquitoes were not blood-fed but were provided with cotton soaked with sucrose solution until 24 h before being released in the experimental huts. Prior to release, five cohorts each with 25 female mosquitoes each were marked with a unique color of fluorescent powder that corresponded to a single, specific experimental hut in which a cohort was released. This marking facilitated monitoring of mosquito movement from release origination to recapture location.

Study design

The study was performed using previously described experimental hut collection protocols [34, 36]. Three trials were performed: i) baseline (no cotton material application); ii) low transfluthrin dosage (0.0025 g/m²); and iii) high transfluthrin dosage (0.005 g/m²). A single trial consisted of five experimental days (i.e. five replicates). Transfluthrin-treated cotton was applied to all interior walls of the center hut (H_T), while solvent-only material (control) was applied to all remaining huts adjacent. Treatments remained fixed throughout a single experimental trial. This treatment allocation was purposive to reflect a semi-urban scenario in a dengue endemic setting whereby an SR product is implemented at only one home in a group of connected houses. Movement during the baseline trial was measured prior to transfluthrin-integrated trials to monitor residual impact of treatment across trials. There was a period of at least two weeks between each trial, in which the huts were thoroughly cleaned and left ventilating with traps removed from windows, doors and eaves.

On each experimental day a single, uniquely marked mosquito test cohort was released inside each of the huts (here named according to the side and distance to the treated hut: H_2L, H_1L, H_T, H_1R, H_2R) at 05:30 hours by five persons. Each individual then remained in the hut, under an untreated bed net, during the 12 h evaluation replicate period to generate host-seeking cues and monitor indoor mosquito knockdown. Knockdown on the floor inside huts was monitored at the top of each hour starting from 06:00 hours until 18:00 hours. Specimens were placed into labeled holding cups, immediately removed from huts and provided access to 10% sugar solution to monitor 24 h mortality. Indoor collectors rotated among huts at the beginning of each replicate day to control for host cue bias.

All window, door and eave exit interception traps were monitored for both alive and knocked down mosquitoes every 30 minutes for a standardized 3-minute collection period per trap. Trap monitoring began at 06:00 hours and continued until 18:00 hours by two-person collection teams (five teams total) positioned outdoors. Collector teams rotated among huts at each sampling period (every 30 min) to control for observer bias and remained at a screened base station positioned 50 m from the experimental huts when not sampling to avoid influencing exit behavior due to outdoor human host cues.

At 18:00 hours, hand-held Prokopack aspirators [37] were used by indoor collectors to recapture remaining mosquitoes inside each hut to calculate loss to follow-up. All recaptured mosquitoes (those from indoor Prokopack aspiration and from exit intervention traps) were held with access to a 10% sugar source to monitor 24 h mortality. Color codes were used to record release origin and location of recapture in a single day.

Model

A continuous-time Markov chain model was developed for the analysis of these data [38]. At any given time, mosquitoes can occupy any one of five huts (transient states: H_2L, H_1L, H_T, H_1R or H_2R) or have experienced one of 15 events represented by the absorbing states: X_2L, X_1L, X_T, X_1R or X_2R for the exit traps in each hut, K_2L, K_1L, K_T, K_1R or K_2R for knockdown in each hut and U_2L, U_1L, U_T, U_1R or U_2R for mosquitoes that were unaccounted for at the end of the experiment and were thus lost to follow-up at some unknown time. The infinitesimal generator matrix A contains the rates at which mosquitoes leave one state to move to another, such that a _ij gives the rate at which a mosquito in state i moves to state j. These rates were assumed to be independent of time or previous trajectories; therefore, the time spent in state i before leaving follows an exponential distribution with mean ${a}_i^{-1}$ with ${a}_i={\sum}_{j=1,j\ne i}^{20}{a}_{ij}$. Note that the rates out of the absorbing states are zero and that, given symmetry in the system, the rates for hut 2L and 2R are equivalent (likewise for 1L and 1R). Subscripts in A indicate the distance from the treatment hut. The 20 × 20 matrix A is defined as

$$ \boldsymbol{A}=\left({\boldsymbol{A}}_H\kern0.5em {\boldsymbol{A}}_X\kern0.5em {\boldsymbol{A}}_K\kern0.5em {\boldsymbol{A}}_U\right), $$

(1)

with

$$ {\boldsymbol{A}}_H=\left(\begin{array}{ccccc}-{q}_2\left(1-{r}_2\right)-{q}_2{r}_2-{k}_2-{u}_2& {q}_2\left(1-{r}_2\right)& 0& 0& 0\\ {}{p}_1{q}_1\left(1-{r}_1\right)& -{q}_1\left(1-{r}_1\right)-{q}_1{r}_1-{k}_1-{u}_1& \left(1-{p}_1\right){q}_1\left(1-{r}_1\right)& 0& 0\\ {}0& 0.5{q}_T\left(1-{r}_T\right)& -{q}_T\left(1-{r}_T\right)-{q}_T{r}_T-{k}_T-{u}_T& 0.5{q}_T\left(1-{r}_T\right)& 0\\ {}0& 0& \left(1-{p}_1\right){q}_1\left(1-{r}_1\right)& -{q}_1\left(1-{r}_1\right)-{q}_1{r}_1-{k}_1-{u}_1& {p}_1{q}_1\left(1-{r}_1\right)\\ {}0& 0& 0& {q}_2\left(1-{r}_2\right)& -{q}_2\left(1-{r}_2\right)-{q}_2{r}_2-{k}_2-{u}_2\\ {}\vdots & \vdots & \vdots & \vdots & \vdots \\ {}0& 0& 0& 0& 0\end{array}\right), $$

$$ {\boldsymbol{A}}_X=\left(\begin{array}{ccccc}{q}_2{r}_2& 0& 0& 0& 0\\ {}0& {q}_1{r}_1& 0& 0& 0\\ {}0& 0& {q}_T{r}_T& 0& 0\\ {}0& 0& 0& {q}_1{r}_1& 0\\ {}0& 0& 0& 0& {q}_2{r}_2\\ {}\vdots & \vdots & \vdots & \vdots & \vdots \\ {}0& 0& 0& 0& 0\end{array}\right), $$

$$ {\boldsymbol{A}}_K=\left(\begin{array}{ccccc}{k}_2& 0& 0& 0& 0\\ {}0& {k}_1& 0& 0& 0\\ {}0& 0& {k}_T& 0& 0\\ {}0& 0& 0& {k}_1& 0\\ {}0& 0& 0& 0& {k}_2\\ {}\vdots & \vdots & \vdots & \vdots & \vdots \\ {}0& 0& 0& 0& 0\end{array}\right), $$

and

$$ {\boldsymbol{A}}_U=\left(\begin{array}{ccccc}{u}_2& 0& 0& 0& 0\\ {}0& {u}_1& 0& 0& 0\\ {}0& 0& {u}_T& 0& 0\\ {}0& 0& 0& {u}_1& 0\\ {}0& 0& 0& 0& {u}_2\\ {}\vdots & \vdots & \vdots & \vdots & \vdots \\ {}0& 0& 0& 0& 0\end{array}\right). $$

The rate q _i signifies the movement rate out of a hut. The direction of this movement depends on r _i (proportion of movement directed to outdoors) and, for H₁, it further depends on repellency p₁ (defined as the proportion of indoor movement from H ₁ , adjacent to H_T, directed away from H_T). The knockdown rate k _i is allowed to vary by hut, whereas the loss to follow-up rate u is assumed to be the same across huts. Hereafter, we refer to the exit rate q_ir_i as x_i (Fig. 3).

Fig. 3

Illustration of experimental hut design and associated model parameters, with q = movement rate, p = proportion of between-hut movement directed away from the treated hut (repellency), r = proportion of movement directed outdoors, q r = x = exit rate, k = knockdown rate and u = loss to follow-up rate. The red hut is the treated hut H_T where the SR treatment is applied. The subscripts indicate whether the parameter applies to H_T (subscript T) or to a hut one or two removed from H_T

The dynamics of the probabilities P _ij (t) of occupying any of the 20 states are governed by a system of differential equations with rates A and are known as the backward Kolmogorov differential equations [38]

$$ \frac{d\boldsymbol{P}}{dt}=\boldsymbol{AP}(t). $$

(2)

From this, we can derive the rates of change in the probability of occupying a given state

$$ {\displaystyle \begin{array}{l}\frac{dH_{2L}}{dt}=\left(-{q}_2\left(1-{r}_2\right)-{r}_2{q}_2-{k}_2-{u}_2\right){H}_{2L}+{p}_1{q}_1\left(1-{r}_1\right){H}_{1L}\\ {}\frac{dH_{1L}}{dt}=\left(-{q}_1\left(1-{r}_1\right)-{r}_1{q}_1-{k}_1-{u}_1\right){H}_{1L}+{q}_2\left(1-{r}_2\right){H}_{2L}+{q}_T\left(1-{r}_T\right)0.5{H}_T\\ {}\frac{dH_T}{dt}=\left(-{q}_T\left(1-{r}_T\right)-{r}_T{q}_T-{k}_T-{u}_T\right){H}_T+{q}_1\left(1-{r}_1\right)\left(1-{p}_1\right){H}_{1L}+{q}_1\left(1-{r}_1\right)\left(1-{p}_1\right){H}_{1R}\\ {}\frac{dH_{1R}}{dt}=\left(-{q}_1\left(1-{r}_1\right)-{r}_1{q}_1-{k}_1-{u}_1\right){H}_{1R}+{q}_2\left(1-{r}_2\right){H}_{2R}+{q}_T\left(1-{r}_T\right)0.5{H}_T\\ {}\frac{dH_{2R}}{dt}=\left(-{q}_2\left(1-{r}_2\right)-{r}_2{q}_2-{k}_2-{u}_2\right){H}_{2R}+{q}_1\left(1-{r}_1\right){p}_1{H}_{1R}\\ {}{\left.\frac{dX_i}{dt}\right|}_{i=2L,...\mathrm{2}R}={r}_i{q}_i{H}_i\\ {}{\left.\frac{dK_i}{dt}\right|}_{i=2L,...\mathrm{2}R}={k}_i{H}_i\\ {}{\left.\frac{dU_i}{dt}\right|}_{i=2L,...\mathrm{2}R}={u}_i{H}_i.\end{array}} $$

(3)

By initializing this system in one hut (e.g. H_2L = 1 and all other states are zero at t = 0), solving for this system of differential equations gives the probability that a mosquito released in a given hut occupies a specific state at time t.

The absorbing states (i.e. X _i , K _i and U _i ) represent competing endpoints in the sense that an individual who enters one of these states is no longer capable of entering any of the other states at some future time. This is reflective of the nature of the experimental data, wherein captured mosquitoes can only be classified as knocked down or exited, but not both. The Markov chain accounts for these competing endpoints vis-à-vis the property that the states are discrete and mutually exclusive. In addition, a mosquito released in 2L can only be knocked down in 2R conditional on having moved there prior to the knockdown event. The absence of non-zero rates to any of the absorbing states from other huts ensures this conditionality.

Likelihoods

To estimate A, we fitted eqn. (3) to the data using a likelihood-based approach. The data collected during the experiments consisted of a set of interval- and right-censored time-to-event data. Outcome measures of interest included exiting (i.e. leaving a space), repellency (i.e. reduced entry into a treated space), knockdown, diversion (defined as the movement to a hut other than the release hut) and loss to follow-up, where exiting, knockdown and loss to follow-up are competing events (i.e. the occurrence of one event precludes the occurrence of another). The cumulative, conditional probabilities for all events observed in the experiment can be directly obtained from the solutions of eqn. (3), as detailed in eqns. (4), (5) and (6).

Interval-censored events

Data pertaining to knockdown and exit events are interval-censored between time points t₁ and t₂, with exit events recorded at 30 minute intervals and knockdown events at hourly intervals. Given model parameter set θ, the probability that a mosquito released in H _rel is observed to be knocked down in hut H at time t₂ is

$$ {\displaystyle \begin{array}{l}\Pr \left({t}_1<T<{t}_2,Y={kd}_H|{H}_{rel},\theta \right)=\frac{F_{kd_H}\left({t}_2|{H}_{rel},\theta \right)-{F}_{kd_H}\left({t}_1|{H}_{rel},\theta \right)}{F_{{\mathit{\operatorname{div}}}_H}\left({t}_2|{H}_{rel},\theta \right)S\left({t}_2|{H}_{rel},\theta \right)}\\ {}\kern10.55999em ={K}_H\left({t}_2|{H}_{rel},\theta \right)-{K}_H\left({t}_1|{H}_{rel},\theta \right),\end{array}} $$

(4)

where F(t) denotes the probability that a specific event (here knockdown and movement to H) occurred in hut H by time t and S(t) denotes the survival function (i.e. the probability that no knockdown, exit or loss to follow-up has occurred by time t). Exit and knockdown events contain indirect information on the diversion event, namely that the mosquito has moved from its release location to the hut where the event took place before the event occurred. This condition, as illustrated by F_div in the denominator of eqn. (4), is implicitly accounted for within eqn. (3); hence, the absence of conditioning in the second part of eqn. (4).

Loss to follow-up

Of mosquitoes that are not retrieved at the end of the experiment, we know that they must have been lost to follow-up at some point between the start and the end of the experiment with probability

$$ \Pr \left({t}_{start}<T<{t}_{end},Y=u|{H}_{rel},\theta \right)=\frac{F_u\left({t}_{end}|{H}_{rel},\theta \right)}{S\left({t}_{end}|{H}_{rel},\theta \right)}=\sum \limits_{i=2L}^{2R}U\left(t|{H}_{rel},\theta \right). $$

(5)

Here, the rate of loss to follow-up u was assumed constant across the huts.

Right-censored data

Mosquitoes retrieved by the end of the experiment are treated as right censored. Namely, the time before knockdown, exit or loss to follow-up would have occurred is longer than the duration of the study, but by how much is uncertain. In addition, we know that the mosquito moved from the release hut to the hut where it was retrieved with probability

$$ \Pr \left(T>{t}_{end}|{H}_{rel,}\theta \right)=S\left({t}_{end}|{H}_{rel},\theta \right){F}_{{\mathit{\operatorname{div}}}_H}\left({t}_{end}|{H}_{rel},\theta \right)={H}_H\left({t}_{end}|{H}_{rel},\theta \right). $$

(6)

Likelihood function

The overall likelihood of the parameters given the data is equal to the product of the probabilities of each individual observation conditional on the parameters. These observations include the number of mosquitoes exited or knocked down during specific time intervals during an experiment for different release huts, event huts and experimental day, as well as numbers recaptured or lost to follow-up at the end of the experiment, resulting in

$$ {\displaystyle \begin{array}{l}L={P}_{multinom}\Big({k}_{exit_{t=i,H=j, rel=k, day=l}},{k}_{kd_{t=i,H=j, rel=k, day=l}},{k}_{hut_{t=750,H=j, rel=k, day=l}},{k}_{u_{t=750, rel=k, day=l}}\mid \\ {}{p}_{exit_{t=i,H=j, rel=k}},{p}_{kd_{t=i,H=j, rel=k}},{p}_{hut_{t=750,H=j, rel=k}},{p}_{u_{t=750, rel=k}}\Big),\end{array}} $$

(7)

with i time points, j event-huts (H_2L, H_2L, H_T, H_1L, H_2R,), k release-huts and l experiment days (1 to 5). Each k_exit, for instance, denotes the number of mosquitoes exiting from hut H observed at time t, by release hut and experiment day. The corresponding probabilities p are derived as detailed in eqns. (4), (5) and (6) and are assumed to be independent of the experiment day.

Model fitting

We used a Bayesian Markov chain Monte Carlo (MCMC) approach for parameter estimation. Using Bayes’ theorem, we define the posterior probability density of the model’s parameters (θ) given the data as

$$ \pi =P\left(\theta | data\right)=\frac{P\left( data|\theta \right)P\left(\theta \right)}{\int P\left( data|\theta \right)P\left(\theta \right) d\theta}, $$

(8)

where P(θ) is the prior probability of the parameters. We utilized beta-distributed priors with median 0.5 for p₁ and median 0.25 for r _i (i.e. a mosquito is three times as likely to move to an adjacent hut than to move outside), a gamma-distributed prior with mean 0.02 on the movement rates q_i (i.e. average time before moving to another hut of 50 minutes) and uniform priors for the remaining parameters (see Table 1 for distribution parameters). Average times before exiting from each hut (1/q _i ) were constrained between 5 minutes and 20 hours and the average time until knockdown (1/k _i ) between 12 hours and 10 days [10, 19]. We explored the parameter space of θ more broadly using the Metropolis-Hastings algorithm.

Table 1 Parameter definitions and prior probability distributions for each

We started from an initial parameter set θ₁, which was randomly sampled from uniform distributions with bounds: q: 360^-1–30^-1, p₁: 0.5–1, r_i: 0–0.5, k: 1400^-1–720^-1 and u: 2000^-1–1000^-1. A new parameter was proposed such that λ₂ = λ₁+Λ, where Λ is a random value from a truncated normal proposal distribution g with mean λ₁ and standard deviation formulated relative to λ₁ and selected so as to ideally have an acceptance rate between 10% and 50% [39]. Which parameter was updated at a given iteration was determined by taking a random draw from a categorical distribution with 11 categories (i.e. the number of model parameters to be estimated) and equal probabilities for each parameter. The probability for λ₂ to be accepted depends on the likelihood of both θ₁ and θ₂ according to the Metropolis-Hastings rule as

$$ A\left({\theta}_1,{\theta}_2\right)=\min \left(1,\frac{\pi_2}{\pi_1}\frac{g\left({\lambda}_1|{\lambda}_2\right)}{g\left({\lambda}_2|{\lambda}_1\right)}\right), $$

(9)

where θ₂ differs from θ₁ only with respect to λ and g denotes truncated normal proposal distributions (between zero and one for each of p and r and from zero to infinity otherwise):

$$ {\displaystyle \begin{array}{l}g\left({\lambda}_2|{\lambda}_1\right)=\frac{P\left(\Lambda ={\lambda}_2\right)}{P\left(\Lambda \le 1\right)-P\left(\Lambda \le 0\right)}\mathrm{for}\ {p}_1\;\mathrm{and}\ {\mathrm{r}}_i\;\\ {}g\left({\lambda}_2|{\lambda}_1\right)=\frac{P\left(\Lambda ={\lambda}_2\right)}{1-P\left(\Lambda \le 0\right)}\ \mathrm{for}\ \mathrm{all}\ \mathrm{other}\ \mathrm{parameters},\end{array}} $$

(10)

where Λ is normally distributed with mean λ₁ and standard deviations corresponding to each parameter’s proposal distribution.

In the event that the acceptance probability was larger than a randomly generated uniform value between zero and one, θ₂ was accepted into the chain. Otherwise, θ₁ was retained. Multiple iterations of this routine were performed (n = 90,000). This process was repeated five times starting from different initial parameter sets to assess convergence using the Gelman-Rubin (GR) statistic [40]. The resulting chains of accepted parameters (φ), after discarding a ‘burn-in’ period (10,000), were combined to represent our sample from the posterior distribution (π).

Simulation experiments

To validate the accuracy of the model-fitting algorithm, we simulated data with a known data-generating process corresponding to our likelihood formulation and with known model parameters. Probabilities for released mosquitoes to occupy a specific state over time were derived using eqn. (3). As follows from eqns. (4–6), these probabilities are defined for interval- and right-censored events. Random draws from a multinomial distribution with the simulated probabilities and a given number of released mosquitoes were taken to simulate numbers of mosquitoes occupying each state at the time points at which sampling was simulated to occur. In general, these simulation experiments were designed to mimic features of the empirical experiments.

Ten distinct simulated parameter sets were used to validate the accuracy of our statistical inference framework. These parameter sets were sampled from across the composite parameter space θ using the Sobol algorithm [41, 42], where the same bounds to this sampling space were applied as for the prior distributions (Table 1). Data were simulated for different numbers of released mosquitoes (25: field scenario; 1000: large sample size scenario) for five replicates per parameter set and fitted to eqn. (3) as described before (n = 60,000, of which 10,000 was ‘burn-in’).

Results

Validating the inference methodology

We first validated the inference framework against data simulated with the system of ordinary differential equations described in eqn. (3), with an observation process that mimicked the field experiment and with parameters reflecting the range of values in the prior distributions.

Large sample size scenario

In the large sample size scenario (five replicates with 1000 released mosquitoes each), we accurately estimated the values of all parameters used in the simulations. All true parameter values fell within the 95% highest density interval (HDI) of the estimated posterior distributions (Fig. 4). Most posterior medians approximated the true parameter well (Pearson’s r > 0.98), but somewhat less so for knockdown in the treated hut (Pearson’s r = 0.74). Posterior distributions were relatively wider for rate parameters associated with the treated hut (x_T and k_T). Standard deviations of these parameters were a fraction (i.e. 11% and 12%) of their respective medians, whereas the s.d.:median ratio was below 3.5% for all other parameters. This reduced precision may be a consequence of the fact that rate parameters associated with huts other than the treated hut were informed by twice as much data as were the rate parameters associated with the treated hut, which derives from our assumption of shared parameters for huts a given distance from the treated hut (Fig. 3). GR statistics were below 1.1 for most simulation sets (average 1.04). When simulation sets resulted in parameters with GR statistics above 1.1, these were related to mosquito movement (q_i, r_i and p₁) and were most commonly associated with the untreated huts (Additional file 1: Table S1). This indicates that those parameters may be among the most difficult to estimate.

Fig. 4

Estimated parameters from simulation experiments for five replicates of 1000 released mosquitoes (large sample size scenarios) with the true value used in the simulation (blue diamonds) and the estimated median (black circles). The dashed gray line depicts p₁ = 0.5, i.e. no repellency effect. Each estimate was based on five chains with distinct starting conditions. 60,000 MCMC iterations were performed inclusive of a ‘burn-in’ period of 10,000

Field scenario

We also tested the performance of the inference framework on data simulated with the same number of mosquitoes as were released in the experiments (Fig. 5). All true parameter values fell within the 95% HDI, but the posterior medians were less consistent with the simulated values (r > 0.8 for all but x_T: 0.68; k₂:0.03; k_T: -0.29) than under the large sample size scenario (Fig. 5). No systematic underestimation or overestimation was observed based on these simulations, suggesting that the additional discrepancy between simulated and inferred parameter values in the field scenario relative to the large sample size scenario was due to stochasticity associated with the smaller sample size in the field scenario (i.e. n = 25 vs n = 1000). GR statistics were, across all parameters and simulation sets, close to 1 (average GR 1.01) (Additional file 1: Table S1).

Fig. 5

Estimated parameters from simulation experiments for five replicates of 25 released mosquitoes (field scenarios) with the true value (blue diamonds) and the estimated median (black circles). The dashed gray line depicts p₁ = 0.5, i.e. no repellency effect. Each estimate is based on five chains with distinct starting conditions. 60,000 MCMC iterations were performed inclusive of a ‘burn-in’ period of 10,000

Product effects on mosquito behavior

We first fitted the Markov chain model to the experimental hut data allowing all parameters to vary. Strong correlations between r _i , q _i and p _i indicated that these parameters were not identifiable given that a wide range of combinations of values of these parameters explained the data equally well (Additional file 1: Figure S1). To resolve this identifiability issue, we fitted the exit rate x _i as a single composite parameter (q _i r _i ). The rate of movement between huts is directly related to the exit rate; namely, it is a proportion (1-r_i) of the overall movement rate out of a specific hut (q _i ). In doing so, we fixed the values of r _i at the medians of the posterior marginal density of the r _i corresponding to each hut that was obtained from the full parameter fit on the baseline data set (r_T: 0.040, r₁:0.036 , r₂: 0.088) (Additional file 1: Figure S1). This reduced the amount of cross-correlation from Pearson’s r as high as 0.86 in the original (Additional file 1: Figure S1) to as low as 0.72 upon fixing r _i (Additional file 1: Figure S2, Figure S3 and Figure S4). Most importantly, it markedly improved convergence from GR statistics as high as 1.38 (q ₂ , low dosage) to as low as 1.00 for all parameters after fixing r _i , indicating that other parameters became identifiable once this adjustment was made (Additional file 1: Figure S7, Figure S8, Figure S9). Choosing either the 2.5th (r_T: 0.031, r₁:0.027, r₂: 0.065) or 97.5th (r_T: 0.054, r₁:0.040 , r₂: 0.11) percentile of r _i instead did not affect this conclusion (Additional file 1: Figure S5 and Figure S6). Acceptance rates for each chain tended to remain relatively constant following a ‘burn-in’ period and varied across chains and parameters within the range of 21–54%.

Exit and movement rates

Under baseline conditions (no chemical), exit rates (x _i ) from huts at different distances i from the treatment hut were relatively similar (medians for x_T: 2.2 ×10^-3, x₁: 1.6 ×10^-3, x₂: 1.8 ×10^-3) (Fig. 6a-c). In subsequent treatment experiments, exit rates out of the treated hut were reduced relative to the baseline in response to both the low (RR = 0.70, HDI = 0.62–1.09) and the high transfluthrin dosage (RR = 0.70, HDI = 0.40–1.06), with no perceptible difference in the respective effects of the two dosages (Fig. 6c). This effect carried over to the adjacent huts (H₁) with exit rates lower than observed in the baseline experiment (low: RR = 0.79, HDI = 0.59–1.01; high: RR = 0.66, HDI = 0.50–0.87) (Fig. 6b). In the huts furthest from the SR application (H₂), the low dosage had no effect on exit rates relative to when no SR was applied (RR = 0.94, HDI = 0.72–1.18). In contrast, the high dosage reduced exit rates (RR = 0.71, HDI = 0.54–0.92) in all huts adjacent to the source of transfluthrin, including the furthest adjoining structures (Fig. 6a). Given that the proportion of movement that was directed outdoors (r _i ) was held constant in this exercise, these results on exit rates (q _i r _i ) are directly proportional to movement rates (q _i ).

Fig. 6

Posterior distributions of model parameters fitted to experimental data for the baseline (gray), low (orange) and high (pink) transfluthrin dosage for the treated hut (subscript T) and huts one or two removed from the treated hut (subscript 1 and 2, respectively). a-c Rates at which mosquitoes exit the huts. d Proportion of movement from H₁ (hut directly adjacent to the treatment hut) away from the SR, where the dashed line indicates p₁ = 0.5, i.e. no repellency effect. e-g Knockdown rates. h Loss to follow-up rates. Under this parameterization, the movement rate q _i is exactly equal to the product x _i / r _i . The algorithm was run for 90,000 iterations inclusive of a ‘burn-in’ period of 10,000

Repellency

In the baseline experiment, mosquitoes moved away from or towards the treated center hut (H_T) with roughly equal probability (p₁ = 0.54, HDI = 0.48–0.59), although with a possible slight preference for movement away from H_T (odds of moving away = 1.16, HDI = 0.92–1.41) (Fig. 6d). In the experiment using low-dosage SR treatment, significant repellency from the treated center hut was observed (odds = 1.64, HDI = 1.30–2.09), with a median probability of moving away from this hut of 0.62 (HDI = 0.57–0.68) (Fig. 6d). In the high-dosage treatment, repellency was still clear (odds = 1.35, HDI = 1.04–1.67), but the effect was somewhat smaller (p₁ = 0.57, HDI = 0.52–0.63) (Fig. 6d).

Knockdown

Knockdown was a very rare event during baseline experiments (2/125 mosquitoes across all five replicates). As a consequence, estimates of knockdown rates in the baseline approached the lower boundary of the prior distribution (medians for H _T = 5.8 × 10^-5_,H₁ = 4.4 × 10^-5, H₂ = 4.0 × 10^-5) (Fig. 6e-g). There was no effect of the low SR dosage on knockdown rates relative to the baseline, both in the treated hut H _T (RR = 1.39, HDI = 0.26–3.84) (Fig. 6g) and in the H₁ huts directly adjacent (RR = 1.00, HDI = 0.45–1.76) (Fig. 6f). In the H₂ huts furthest away from the treatment, a somewhat increased knockdown rate was observed in response to the low dosage relative to the baseline (RR = 1.37, HDI = 0.64–2.46) (Fig. 6e). Knockdown rates in the high-dosage scenario were elevated in all huts, in particular in the H _T treatment huts (RR = 8.37, HDI = 2.11–17.35) (Fig. 6g) but also in the H₁ and H₂ huts (H₁: RR = 1.39, HDI = 0.52–2.69; H₂: RR = 2.22, HDI = 0.96–3.86) (Fig. 6e, f).

Loss to follow-up

Rates of loss to follow-up were similar across the baseline and two SR treatment experiments, although there was a signal for a small increase in these rates with increasing dosage (low: 5%; high: 8%) (Fig. 6h). In comparing posterior samples across dosages, a signal for a positive dose-response relationship [i.e. u(high) > u(low) > u(baseline)] was confirmed in 61% of samples from the posterior. Loss to follow-up rates in the baseline scenario were an order of magnitude higher than the estimated knock-down rates. While some of mosquitoes that were lost to follow-up at the end of the experiment could have been knocked down, the differences in trends between treatment effects on knockdown and loss to follow-up indicate that this could likely only be true for a small portion of the mosquitoes.

Time spent in a hut

The total amount of time a mosquito spent in each hut results from the composite of treatment effects. By running simulations of the system of ordinary differential equations (eqn. (3)) with the estimated posterior parameter values, we derived a posterior estimate of the proportion of the time a mosquito spent in each hut relative to the total time a mosquito was in the hut system (i.e. before exit, knockdown or loss to follow-up). This proportion was found to be similar but slightly reduced for the treated hut H_T relative to the baseline scenario in either treatment scenario (Fig. 7e) and without any effect in the downstream huts H₁ and H₂ (Fig. 7a, c). However, when considering the total duration of the experimental day, the proportion of time spent in the adjacent, downstream huts H₁ and H₂ was higher during experiments using both low and high SR dosage than during baseline (Fig. 7b, d). This was a result of reduced exit rates and thus an overall increase in time spent in the hut system as a whole (Fig. 6a, c).

Fig. 7

Distributions of time mosquitoes spent in each hut relative to the posterior median of the baseline scenario (gray), low (orange) and high (pink) transfluthrin dosage for huts two (a and b) or one (c and d) removed from the treated hut (e and f). The left column signifies the proportion of time spent in each hut before having experienced an event (a, c and e), where kd is knockdown. The right column signifies the proportion of the total experiment time spent in each hut relative to the baseline (b, d and f)

Discussion

Novel Ae. aegypti vector control strategies are currently being evaluated to address challenges related to dengue transmission expansion [2]. Spatial repellent (SR) products, which release volatile chemicals into treated spaces to interrupt host-vector contact, are among these [14]. One challenge for evaluating the efficacy of SRs, and other products that may include non-lethal outcomes, is characterizing the multifaceted, incapacitating effects of a given product on mosquito behavior under field conditions and across a dosage gradient, which is critical for identifying opportunities to exploit new mechanisms of action and guiding product formulation and implementation under operational conditions. To quantify concurrent and downstream (i.e. untreated spaces adjacent to the treated space) effects of a transfluthrin-based SR product on Ae. aegypti behavior, we used a continuous-time Markov chain model informed by experimental data under a Bayesian inference framework.

Examination of posterior estimates of model parameters showed that test mosquitoes were deterred from entering the experimental hut where the product was located and that this effect was stronger for the lower dosage SR, an effect of transfluthrin consistent with other studies [43, 44]. Posterior estimates of model parameters also indicated lower movement rates out of the treatment hut (either to a neighboring hut or out of the hut system) under both treatments, presumably due to confusion effects, similar to previously demonstrated effects of metofluthrin [13]. Under the higher SR dosage, the reduced out-of-hut movement was noticeable in adjacent untreated huts as far as two huts away from the SR application. This dose-dependent effect may be a result of irritancy and disorientation, as reviewed in Bibbs & Kaufman [45]. Similarly strong effects in adjacent connected huts were observed on knockdown rates, which were markedly increased in all huts in the presence of the higher SR dosage, yet little effect on knockdown was observed at the lower dosage. These results are in agreement with dose- and distance-dependent effects on knockdown observed upon metofluthrin exposure in Ae. aegypti [46].

While metrics of repellency and exit rates from traditional experimental hut studies are often affected by treatment effects on mosquito knockdown, our inference framework explicitly accounts for such confounding effects. We found knockdown and loss to follow-up rates to be small relative to mosquito movement rates, indicating that for the formulations considered, effects on repellency and exit rates could have meaningful impacts on mosquito-human contact rates.

We validated our inference method by demonstrating its ability to accurately estimate the model’s parameters given simulated data. This assessment was conditional, however, on the assumption that the model is a realistic representation of reality. Some of the known limiting assumptions of our analysis include (i) effects that depend on distance from the treated hut rather than on each hut individually; (ii) equal loss to follow-up across huts; and (iii) time-invariant parameters. Of these, the first may be most problematic when considering that air flow within the hut system could result in asymmetric effects of transfluthrin dispersion to huts of the same distance from the treated hut but on different sides of it [47]. In principle, it would be possible to account for such factors in future studies by measuring air flow and incorporating its effect on the data through an appropriate modification of the model. For example, repellency (p ₁ ) could be allowed to vary across huts and treated as a function of readings from a wind gauge. Posterior estimates of the parameters governing the relationship between wind and p ₁ would then allow for inferences about the repellency of the product under varying airflow conditions and beyond those observed in the current experiment. Planning for required sample sizes and sampling schemes for such experiments would benefit from our model and results by using our posterior distributions to inform prior distributions in those future studies [48].

Repellency (i.e. reduced entry into the treated hut from adjacent connected huts) and increased knockdown reduced the overall time that mosquitoes spent transiently in the treated hut, whereas decreased movement rates have the potential to offset this effect. The result of the SR’s impact on the time mosquitoes spent transiently in the treated hut indicates potential for such a product to limit human-vector contact (and thereby reduce the probability of pathogen transmission) in the treated hut. Furthermore, the reduced exit rates that we observed may be a result of confusion and disorientation and may be accompanied by impaired host-seeking [45]. However, because blood-feeding was prohibited in our study, it is uncertain to what extent host-seeking and blood-feeding behaviors of these mosquitoes exposed to the SR may have been affected in the current study. Other studies using similar volatile products have shown these effects to also be associated with reduced human landing [13, 49, 50]. The inclusion of blood-feeding metrics in experiments with volatile pyrethroid products using anophelines under field conditions [51] and against the topical repellent DEET using Ae. aegypti in the laboratory [52] have been valuable in establishing expectations of such synergistic chemical effects.

The effect of SR products on untreated neighboring premises has been a consistent and critical question to the public health value of these products [14, 53]. Three aspects of our results suggest that the risk of diversion (i.e. movement of mosquitoes from a treated space to an adjacent untreated space) may be limited for the formulation used in our experiments. First, SR exposure reduced movement rates between huts. Secondly, there was a marked increase in knockdown in untreated huts at the high SR dosage. At the same time, there was also a marked reduction in exit rates out of untreated huts, which resulted in prolonged time spent in adjacent huts. Evaluating the overall potential for diversion based on these effects will require pairing experimental results such as ours with theory that is capable of translating this range of behavioral effects into estimates of their epidemiological consequences [15, 54, 55].

Under our experimental design, we cannot distinguish between downstream effects caused by volatile particles dispersed into untreated huts or by a residual, post-exposure effect of transfluthrin on mosquitoes that are exposed in the treated hut and move elsewhere. Indeed, this highlights that SR post-exposure effects remain understudied. A mosquito that remains alive in a treated space, and/or does not exit after being exposed to a SR product, may exhibit attraction-inhibition to human hosts inside the space or may be inhibited to blood-feed due to being chemically incapacitated. Underlying physiological mechanisms causing these responses are yet to be characterized. Correlations between air sampling measurements in experimental huts and mosquito behavior responses have been explored in previous studies using spatial repellents but with limited success due to limits of chemical detection and quantification [56,57,58]. Combining air chemistry inferences of specific active ingredients (i.e. vapor pressure or particle weight) with environmental data (i.e. air current or flow rate) into our new inference framework is therefore warranted and could enable quantification of the extent to which downstream effects result from movement of the volatile chemical or movement of exposed mosquitoes with lingering post-exposure effects. The latter possibility has been indicated in other studies to have potential for innovative applications of SRs [59].

Conclusions

The need for development and efficient testing of new vector control products and innovative formulations of existing tools is evident [2, 60]. Advancing the characterization of SRs and other volatile vector control products must include understanding how all product-induced effects, independent or in combination, contribute to potential impact on pathogen transmission. Here, we advance SR product characterization methods by integrating multiple aspects of the complexity of assumed product-induced responses. The complementary experimental design and inferential framework we introduce provides a flexible approach for estimating a product’s effects on mosquito behavior in a quantitative and probabilistic fashion. Using this framework, we disentangled complex, dose-dependent effects of transfluthrin on mosquito behavior. While a low-dosage treatment resulted in notable repellency and minor knockdown in Ae. aegypti mosquitoes, higher dosages presented with much higher knockdown rates yet minor repellency effects. Exposure to transfluthrin has the potential to diminish mosquito movement, host attraction-inhibition and/or blood-feeding inhibition, possibly due to confusion or disorientation of the mosquito. The use of an experimental row-house hut design allowed for the estimation of transfluthrin effects beyond a treated space, which is of particular interest as adjacent, inter-connected houses are common in semi-urban environments where SR products are under evaluation as a tool against arbovirus transmission. These findings highlight the need to estimate context and dosage-specific diversion (product-induced movement to an untreated space) and downstream (adjacent spaces) effects under scenarios similar to operational settings for which the product is intended to be used.