A Spatial Model of Mosquito Host-Seeking Behavior

Mosquito host-seeking behavior and heterogeneity in host distribution are important factors in predicting the transmission dynamics of mosquito-borne infections such as dengue fever, malaria, chikungunya, and West Nile virus. We develop and analyze a new mathematical model to describe the effect of spatial heterogeneity on the contact rate between mosquito vectors and hosts. The model includes odor plumes generated by spatially distributed hosts, wind velocity, and mosquito behavior based on both the prevailing wind and the odor plume. On a spatial scale of meters and a time scale of minutes, we compare the effectiveness of different plume-finding and plume-tracking strategies that mosquitoes could use to locate a host. The results show that two different models of chemotaxis are capable of producing comparable results given appropriate parameter choices and that host finding is optimized by a strategy of flying across the wind until the odor plume is intercepted. We also assess the impact of changing the level of host aggregation on mosquito host-finding success near the end of the host-seeking flight. When clusters of hosts are more tightly associated on smaller patches, the odor plume is narrower and the biting rate per host is decreased. For two host groups of unequal number but equal spatial density, the biting rate per host is lower in the group with more individuals, indicative of an attack abatement effect of host aggregation. We discuss how this approach could assist parameter choices in compartmental models that do not explicitly model the spatial arrangement of individuals and how the model could address larger spatial scales and other probability models for mosquito behavior, such as Lévy distributions.


Introduction
The transmission of infectious agents is heterogeneous in the sense that the risk for infection and infectiousness are unevenly distributed over the population [1]. This heterogeneity is an important factor in predicting the spread of an infection; ignoring it may result in misleading inference about transmission dynamics by underestimating the probability that an infectious agent will persist [2]. Sources of heterogeneity in mosquito-borne disease transmission include uneven distribution of hosts and breeding sites [3], differences in host susceptibility to infection [4], and the nonuniform distribution of mosquito bites among spatially distributed hosts [5] (per-capita biting rates). The number of mosquitoes that bite a host depends, among other things, on the number that find it. Host finding by mosquitoes is largely driven by olfactory cues that are given off by individual hosts [6]. The spatial arrangement of hosts is likely to affect the spatial distribution of the odor plume and thus the mosquitoes' ability to locate and feed on them. For example, it may be easier for a mosquito to find a large roost of birds than to find an single bird. Unless the probability of finding a host is exactly proportional to the density of hosts, unevenly distributed contact rates on individual hosts, and thus heterogeneous disease transmission, will result. Understanding the dynamics of odor-driven mosquito-host interaction is fundamental to a detailed mechanistic understanding of mosquito-borne transmission.
Despite its important epidemiological implications, experimental data on the relationship between host aggregation and mosquito-host contact rates are sparse, and thus, little is known definitively about mosquito strategies for finding hosts. The goal of this work is to propose a modeling framework that can provide preliminary insights into the relative effectiveness of different hostseeking strategies used by mosquitoes. In this paper, we develop, assess, and utilize a mathematical model to simulate both the odor plume in the presence of wind and the host-seeking behavior of mosquitoes in response to that odor plume. Mosquitoes are modeled as discrete agents that fly continuously in search of discrete hosts (birds). The flight direction and speed of individual mosquitoes is influenced by wind and odors emitted by the hosts. The wind is composed of a deterministic, large-scale component plus a stochastic component to represent small-scale eddies and fluctuations. In this study we focus on mosquito host-seeking behavior within meters of a host, motivated by indoor experiments involving mosquito-bird interactions [7]. The Discussion section describes the potential for larger scale simulations that incorporate more of the host-seeking process. We hope that the model may serve as a virtual laboratory for testing different hypotheses for how mosquitoes use odor plumes to locate potential hosts and may offer guidance both to the planning of experimental studies and to the interpretation of experimental data.
This contribution joins a group of other important mathematical models for simulating the host-seeking dynamics of mosqui-toes. Commonly, models are based on parameters such as contact rates and probabilities (of feeding, diversion, dying, etc.). Often these models assume homogeneous mixing and do not include space explicitly, so that the probabilities and contact rates are independent of location. In [8], the authors describe a kinetic, non-spatially explicit model to determine how much coverage is enough to protect individuals who do not use insecticide-treated nets (ITNs) in the context of malaria transmission. The model has also been modified to explore the effects of bed nets [9] and to introduce 'bloodless' hosts such as baited traps [10]. The latter also modeled the host-seeking process as a non-host oriented kinesis followed by a host-oriented taxis once the mosquito has encountered an odor cue. A related model for the mosquito feeding cycle addressing the effect of ITNs on malaria transmission is found in [11]. One model that does include space explicitly is that in [12] for plume finding in the context of an underwater substrate. The authors use a model in which the probability of detecting the plume depends on the seeker's movements, on the fluid motion, and on the plume shape in space and time. They report that as plume-detection success increases, the efficiency decreases (the process is slow). Our results are consistent with this observation.

Mosquito host-seeking behavior
We provide a review of the biological literature that influenced our modeling choices. The scenario that we consider is illustrated in Fig. 1, where hosts are distributed in groups, or patches, and they emit an odor (e.g. CO 2 ) that gets carried by and diffused in the wind in the same way that a puff of smoke dissipates in time. Uniform, laminar wind extends host odor into a long, thin plume with sharp transverse gradients and shallow longitudinal gradients. If the wind is turbulent, the odor plume is highly intermittent, but still retains relatively shallow average longitudinal gradients compared to the transverse gradients [13]. We class mosquito host-seeking behavior as either plume finding, which is flight in search of an odor plume, or plume tracking, which is flight within the odor plume (these terms are adopted from [12]).
Plume finding: Absence of odor cues. There is little consensus in the literature about mosquito host-seeking behavior in the absence of odor cues. If wind is absent, orientation may be determined by large visual features in the environment [14] or may be characterized as a directionally unbiased random walk, also called kinesis [10]. If wind is present, then mosquitoes may deliberately choose to fly upwind, downwind, or crosswind in search of a host, which are types of directional searching referred to as anemotaxis [15].
Each of the three anemotactic behaviors has been described as plausible based on either experimental or theoretical work. Mosquitoes typically fly upwind in laboratory wind tunnels even when there is no odor present [16,17]. In wind tunnel experiments, a low velocity artificial wind is blown across an odor source down an enclosure toward a mosquito entrance. Mosquitoes are released into the tunnel and their flight path is videotaped.
In control experiments where no odor is added to the wind, many mosquitoes fly upwind and some even locate the ''source'' -the wind entrance into the tunnel. However, a set of field experiments reported in [18] provided evidence for downwind flights in hostseeking Mansonia spp. In these experiments, a human subject was surrounded on the downwind side by a tall, hemispherical fence of radius 18 m, to exclude mosquitoes using an upwind search strategy. By comparison with controls, the authors concluded that some mosquitoes employ a downwind plume-finding strategy. In [19] it is argued using geometry that crosswind searching is the most effective when the odor plumes are long and thin. If the variability of the wind direction is greater than 30 degrees, then a mathematical argument shows that upwind or downwind searching is optimal [20].
Mosquito response to wind depends on the strength of the wind. A typical mosquito flight speed is 1 m/s [15,17]. As reviewed in [21], mosquitoes fly faster than the wind speed when they are in the low velocity boundary layer that forms adjacent to the ground. If they ascend too far (especially during the daytime when turbulence is greater), they are swept up and transported passively for long distances. Some mosquitoes may have adapted to use this as a deliberate migration mechanism. However, most of the flights that mosquitoes make are short-range, appetitive flights seeking a blood meal or oviposition site near the mosquito's home territory. Appetitive flights are disrupted by sufficiently high wind speeds, although the wind speeds at which this happens vary by species and geographic region. Wind speeds as low as 0.8 m/s (3 km/h) have been reported to drastically reduce the number of mosquito host-seeking flights, but no reduction in mosquito flight was documented in other situations for speeds as high as 3-8 m/s (11-29 km/h).
Plume tracking: Presence of odor cues. When a mosquito encounters an odor plume, it uses the odor plume and the wind to guide its flight to locate the host. Odor cues are complex olfactory signals released from a host's skin and breath. Many compounds are known to excite the chemoreceptors of mosquitoes. CO 2 , for example, is an important component in the odor plume that activates and helps maintain plume tracking [17,22,23]. Laboratory wind tunnel experiments indicate that sustained flight only occurs in the presence of an intermittent CO 2 signal and not in uniform concentrations [23]. Other wind tunnel experiments suggested that broad, well-mixed CO 2 plumes may inhibit upwind flight while turbulent plumes of the same concentration induce upwind flight [16,24,25]. The importance of lactic acid, various aldehydes, and whole host odors for the location of hosts by mosquitoes has also been demonstrated [16,22,25,26]. Dekker et al. [25] note that there is not a single model that fits mosquito flight response to all relevant host odors, even within a single species.
Mosquitoes exhibit different behavior in windy and windless conditions. In the absence of wind within the odor plume, mosquitoes must rely solely on odor cues [13] or on large features in the visual environment [14]. Mosquitoes probably estimate the direction of odor increase (the gradient) by the mechanism of klinotaxis, as is conjectured for tsetse flies [27]. During klinotaxis, an

Author Summary
Mosquito-borne diseases can spread when a mosquito bites a vertebrate host to obtain a blood meal for egglaying. The first step in the transmission process consists of the mosquitoes seeking and finding a host. Mosquitoes use the wind direction and odors, such as carbon dioxide, emitted by the hosts in order to locate a host to bite. We present a spatial computational model of the host-seeking process in a region where hosts are heterogeneously distributed in clusters. The model is used to analyze the success in finding hosts once the mosquitoes are close to the host. We show that the number of mosquito-host contacts increases as hosts become more widely spaced within their clusters; that mosquito flight perpendicular to the wind leads to greater success in locating a host; and that the number of bites per host decreases when hosts aggregate into larger clusters.
organism samples the host odor at one location, then moves and samples it again, using its memory of the concentration to choose the next direction [13]. Laboratory and field wind tunnel experiments suggest that under windy conditions within the odor plume, the mosquitoes travel upwind to locate the source [16,24,28]. Mosquitoes appear to infer wind direction from the optical flow of ground features relative to their position [27]. It is known that many organisms have characteristic turns in their upwind flight path (e.g. moths [13,27]), but mosquitoes exhibit highly irregular upwind flight [29].

Model
We formulate our model for the intermediate range flights of night-active mosquitoes such as Culex quinquefasciatus feeding on roosting birds, incorporating the biological information from the previous section. We construct a model that accommodates different plume-finding and plume-tracking mosquito behaviors, wind velocities and host locations. The model is two-dimensional and all of the dynamics are assumed to take place at a fixed distance from the ground.

Odor plumes
In the model, hosts emit a single gaseous compound that attracts mosquitoes, is convected by the wind, and diffuses in the air. For the purpose of this paper, we assume that the attractant is CO 2 . The CO 2 distribution over time is modeled by a convectiondiffusion partial differential equation. The convection velocity of the wind is given by a vectorṼ V (x,y,t), and the concentration C(x,y,t) of CO 2 is described by the equation where t is time and (x,y) are spatial coordinates on a square domain of length L. Throughout the paper, we refer to this square computational domain as the ''CO 2 simulation region.'' The constant diffusion coefficient D reflects the rate at which CO 2 Figure 1. Computational domain. This is a schematic of the CO 2 simulation region where the odor plume is computed. At the edges of the simulation region, the CO 2 concentration is carried out by the wind. The patches represent smaller subregions where the hosts are located. The spatial distribution of hosts within each patch is uniform and there can be multiple patches. The hosts are stationary in the following simulations, although this is not a limitation of the model. The mosquitoes are initially placed in a subregion inside the CO 2 simulation region but are allowed to leave it (and possibly reenter it) during the simulation. The wind exists everywhere, even outside the CO 2 simulation region. doi:10.1371/journal.pcbi.1002500.g001 diffuses in air in the absence of wind. The term C s (x,y) represents the odor cue emitted at a constant rate by the hosts in units of ppm/s: C s (x,y)~J 0 , at host locations 0, elsewhere: The wind velocityṼ V consists of two components: whereŨ U is used to introduce drifts or relatively large features produced by the air andŨ U r is a stochastic velocity vector used to approximate the effect of small-scale wind variations in the domain. The direction ofŨ U r at each point in the domain is chosen uniformly from ½0,2p) and the magnitude is chosen from a normal distribution centered around zero. The large-scale velocity field U U~½0,U 2 is a constant speed flow from bottom to top (see Fig. 2) in most of our simulations. However, we sometimes use a meandering plumeŨ U m in place of the straight plumeŨ U, which is given by the expressioñ where the frequency a is p=2. It is easily checked that this flow is incompressible.
Equation (1) is numerically evolved using second-order centered differences to approximate the Laplacian and a first-order conservative upwind finite difference method for the convection term (similar to page 636 of [30]). The normal components of the concentration gradient are zero at the boundary (Neumann conditions). We integrated the equations with a forward Euler method and confirmed that our solution converged by verifying that significantly decreasing the spatial grid size and time step size had a minimal change in the solution.
In most of the simulations in this paper, we consider length scales and time frames consistent with the mosquitoes being in close proximity to the hosts. The length of a side of the square domain L is 10 m and the simulations cover time periods of 50-500 seconds. The hosts are situated in the middle of the domain, 5 m from the top and bottom domain edges (see Fig. 2). Initially, the domain is bare of CO 2 and it takes about 45 seconds for the plume to reach the domain edge. At that point, the mosquitoes are released into the domain, with starting positions dependent on their particular flight behavior. For upwind plume-finding behavior, the mosquitoes are all released downwind of the hosts; for downwind and crosswind behaviors, the entry is along the upwind side of the domain.

Mosquito behavioral rules
Mosquitoes and hosts are modeled as discrete individuals, or agents. Hosts are stationary, motivated by the interaction between nocturnal Cx. quinquefasciatus mosquitoes and their roosting bird prey. In this section we explain the mosquito navigation model, which differs during plume finding and plume tracking. Assumptions about mosquito agent behavior include: N The mosquitoes do not affect each other. N The mosquito motion is not restricted to the CO 2 simulation region and may leave and re-enter the region during their random walks.
N Below a CO 2 threshold C 0 , the mosquitoes navigate only using wind direction and move according to a predetermined upwind, downwind, or crosswind pattern. This is plumefinding behavior. N For concentrations above the threshold C 0 , the mosquitoes respond by changing to plume-tracking behavior: moving upwind and toward larger levels of CO 2 . There is a saturation concentration level C sat above which no further changes in concentration can be detected. N When a mosquito comes within a predetermined radius r c of a stationary host, the mosquito is removed from the simulation and a ''contact'' is recorded. In this context, a contact means an attack on the host, regardless of whether it results in a blood meal, a diversion, or death [10]. The term ''contact rate'' refers to the number of contacts per time period, usually the length of one simulation.
N Mosquito agents may move toward higher concentration levels by one of two mechanisms -temporally sampling the concentration level or directly sensing the spatial gradient of the concentration (klinotaxis and tropotaxis respectively in chemotaxis literature [13]).
Mosquito flight direction based on CO 2 concentration (klinotaxis). During plume-tracking behavior, the mosquito chooses a flight segment direction by comparing the current CO 2 concentration level with the concentration previously encountered (CO 2 levels are interpolated to the mosquito location). If the CO 2 level is higher, the mosquito is biased to continue in the same direction. If it is lower, the mosquito is biased to turn around. So the target direction of the mosquito, h, is either plus or minus its previous direction. Given that there is inaccuracy in the ability for a mosquito to sense CO 2 levels, we add stochasticity to the ultimate choice of direction by constructing a window of width 2a centered around h, ½h{a,hza (see Fig. 3).
The quantity a is dependent on the magnitude of the CO 2 change. To calculate it, we set a maximum concentration change that can be sensed, calling it DC sat , and a minimum concentration change DC 0~b0 DC sat , with 0vb 0 v1. Concentration changes below DC 0 are imperceptible. We define a relative concentration change by b~DC(x n ,y n ){C(x n{1 ,y n{1 )D=DC sat , where (x n ,y n ) is the mosquito location at time t n . Then a is computed from where a max and a min are the maximum and minimum allowable window half-widths and F is a ramp function that takes values between 0 and 1: The mosquito direction influenced by the CO 2 concentration, h c , is chosen randomly from a uniform distribution in the window interval: The interval size is largest when the sensory input is lowest, leading to greater uncertainty in direction choice. We refer to this rule as the ''sampling method'' in the remainder of the paper.
Mosquito flight direction based on CO 2 gradient (tropotaxis). Gradient chemotaxis models use the concentration gradient, +C~(LC=Lx,LC=Ly), as the variable sensed by mosquitoes or other chemotactic agents [31,32]. This model is incorporated into our framework by computing the CO 2 gradient on the grid and interpolating it to the mosquito location. The direction of the gradient is the target direction h for the mosquito's next flight segment. Following the procedure in the previous section, we set G sat as the maximum and G 0~b0 G sat as the minimum gradients that can be sensed, and define the relative gradient at the mosquito location by b~D+CD=G sat . Mosquito direction is then computed as in Eqs. (3)- (4). We refer to this as the ''gradient method.'' Mosquito flight direction based on wind velocity. Mosquito agents have a preassigned plume-finding strategy: downwind, upwind, or crosswind. The target angle for wind, denoted h v , depends on the wind direction at the mosquito location, which is evaluated from the large scale windŨ U plus a random component interpolated from the grid in the CO 2 simulation region. For example, if wind is blowing from South to North, a mosquito with a downwind strategy would have a target angle pointing North; a mosquito with an upwind strategy would have a target angle pointing South; and one with a crosswind strategy would have a target angle pointing either East or West (Fig. 2).
The size of the angle window is chosen according to the strength of the windṼ V sensed by the mosquito. We assume that mosquitoes can distinguish wind speeds up to a saturation value V sat . The direction is chosen randomly from a precision window where . The threshold speed for mosquito response to wind is a percentage of the saturation value, V 0~b0 V sat . During crosswind plume-finding behavior, the duration of travel in one direction is a parameter in the model. We choose the duration to be uniformly selected from an interval T cwd~½ T min ,T max . Mosquito flight speed. A mosquito agent chooses a flight speed s between a minimum S min and a maximum S max depending on the local CO 2 level. Unlike the random choice of direction, the speed of the mosquito is defined by a deterministic formula similar to Eq. (3): where the variable b represents the scaled CO 2 concentration C=C sat at the mosquito location. If the concentration is below the threshold C 0 , then the mosquito speed is S max . Similarly, if the concentration is above C sat , then the segment speed is S min . This allows mosquitoes to remain near areas of high CO 2 concentration. A change in speed due to the presence of a chemical signal is properly called orthokinesis [33].
Mosquito position updating. The direction of the flight segment and the speed of each mosquito are updated every DT time units. The updated mosquito position is calculated by where (x n ,y n ) is the mosquito position at time step n. The last term is passive convection of the mosquito in the ambient wind. The quantityd d is a direction vector that varies between plume finding and plume tracking: d d~c where h w from Eq. (5) is based on the wind and h c from Eq. (4) is based on the CO 2 . In plume finding, when there are no odor cues, c~1, and the flight direction depends only on the assigned plumefinding behavior of the mosquito. During plume tracking, c~1=2 and the flight direction is the average of the directions chosen from the concentration and from an upwind strategy. Mosquitoes that have a downwind or crosswind plume-finding behavior will begin to fly upwind in the presence of CO 2 .

Assessment of the model
We assess model performance by comparing different mosquito navigation strategies and by evaluating the sensitivity of the model to a subset of the simulation parameters. We find that the gradient and sampling methods can exhibit comparable performance and that, in general, the parameter set in Table 1 is robust to small changes.
Comparison of simulated host-seeking mechanisms. One might expect that perfect knowledge of the gradient would result in better performance than perfect knowledge of the CO 2 concentration. However, the model includes a window around the gradient direction from which the direction of motion is chosen randomly; in other words, the knowledge of the gradient is imperfect. We find that sufficiently imperfect knowledge of the gradient direction can result in overall behavior that is comparable to the sampling method (with different parameter values). This is an important observation since many existing models use the popular Keller-Segel approach that assumes knowledge of the concentration gradient [31,32].
We placed nine hosts in a regular square grid at 1 ft (0.3 m) intervals in the center of the simulation region in the presence of a meandering velocity field. We examined upwind, downwind, and crosswind plume-finding behaviors, each in combination with the gradient and sampling methods for CO 2 sensing. For the gradient method, we used the parameters a min~p =6, G 0~0 :001 (40 ppm/ meter), and G sat~0 :198 (7900 ppm/meter). The last two parameters are analogous to DC 0 and DC sat for the sampling method. See Eq. (2) for the formula governing the wind and Table 1 for all other simulation parameters. For convenience, all numerical simulations were computed using dimensionless variables, but we report parameters in dimensional units.
We also simulated a random walk strategy (with equal probability of stepping in all directions) in which the odor plume is not sensed at all and the wind has no effect. This mimics mosquito dispersal (mathematical diffusion) independent of the chemical concentration and wind. The comparison of the chemotactic rules to this random walk provides a baseline for determining how different two rules are. The random walk mosquitoes started at the same location as the downwind strategy, but chose a direction at every time step from a uniform distribution on ½0,2p). Their speed was constantly S max , the same speed as the plume-finding behavior of the host-seeking mosquitoes. There are many other possible choices for random walks; see [12].  The proportion (P) of mosquitoes that found a host after 150 s is shown in Fig. 4 for all the mosquito heuristics discussed above. The values shown in the bar graph are averages over 15 simulations in which the time series of random velocity fields is fixed, but the mosquito behavior stochastically varies. The thin error bars represent +2 standard deviations. The random walk (black bars) had by far the fewest contacts, while the gradient and sampling methods have overlapping error bars. After 500 s, there was almost no change in either of the host-seeking methods, but the proportion of mosquitoes that found a host during the random walk increased to 0.25.
For these particular parameter choices, the gradient method and sampling method show similar results and both are substantially different than a random walk. There are enough free parameters in our model of host-seeking so that behavioral heuristics based on either spatial gradient sensing or temporal sampling can deliver approximately the same results. Since mosquitoes almost certainly navigate using klinotaxis [27], we will use the sampling method in the remainder of this work. However, we emphasize that our observations support the use of the Keller-Segel model of gradient sensing as in [31,32].
Sensitivity analysis. The goal of this paper is to characterize the effect of host aggregation and mosquito plume-finding behavior on the mosquito-host contact rate for the parameters in Table 1. In this section, we consider a single host distribution and locally vary a subset of our model parameters independently to assess the robustness of the model with respect to our target parameter set. We find (with a few exceptions) that the parameter set we use with our model is robust to small changes in input parameters.
We simulated crosswind, downwind, and upwind plume-finding behaviors in the presence of a straight odor plume produced by an arrangement of nine hosts spaced 1 ft (0.3 m) in a square patch in the center of the domain. For each plume-finding behavior, we varied each of the ten starred parameters in Table 1 by +10% while holding all other parameters constant. To capture the average mosquito behavior, the simulations were repeated 15-20 times for each plume-finding behavior (upwind, downwind, and crosswind), until a mild convergence criterion was satisfied. In each set of simulations, we used the same sequence of random velocity fields and the same set of random numbers for mosquito behavior, so that the differences between sets were due solely to the parameter perturbations and not to the stochastic nature of the agents or wind. We analyzed two output variables: N the proportion of mosquitoes that find a host anywhere in the domain, P, and N the average time to locate a host, T avg .
When calculating T avg , we exclude mosquitoes that never locate a host.
A local sensitivity index, SI, is a partial derivative of an output variable with respect to an input parameter that is scaled to allow comparisons across variables. High sensitivities are synonymous with large rates of change, while low sensitivities (and model robustness) are associated with small derivatives. The sign of SI indicates if the change is in the same direction (positive) or the opposite direction (negative) as the perturbation. To calculate SI values, we used a centered difference approximation to the partial derivative multiplied by the baseline value of the input parameter over the baseline value of the output variable. Symbolically, if we let I be an input parameter, we have where P { ,P, and P z are the values of the output variable corresponding to the varying input parameter: I{dI, I, and IzdI, respectively. Since we consider variations of +10%, we choose d~0:1. We estimated the error in this method by taking half of the difference between the one-sided partial derivatives: We show the combinations of plume-finding behavior, input parameter, and output variable with the highest sensitivity indices in Table 2. SIs that are small compared to 1 in absolute value indicate a robust response to parameter variation. In the first row, SI~{1:8 means that if the maximum speed of a mosquito (S max ) engaging in upwind plume finding increases by 10%, then the average time for a mosquito to find a host will be decreased by an estimated 18% with error bounds of +3%. This is a large response and represents high sensitivity.
In Table 2, there are two SIs very close to 0.3, and five greater than 0.6. The values greater than 0.6 represent moderate to high sensitivities, while the two values near 0.3 are reasonably low. All combinations that are not shown in Table 2 have DSIDv0:3. The total number of combinations tested were 20 per plume-finding behavior, or 60 total. We see low values of SI in most of the combinations tested, indicating that our model is locally robust to the parameter set in Table 1, although we do not characterize variability due to stochastic effects. The biggest exception to this robustness is the maximum flying speed of the mosquitoes (S max ), which strongly affects T avg in upwind and downwind behaviors and moderately affects P in crosswind plume-finding behavior. Overall, the average time to locate a host, T avg , is more sensitive to input changes than the proportion of mosquitoes that find a host, P.

Results
We use the model introduced in the previous section to address several questions pertaining to space-dependent small-scale plume encounter and host localization. In the following subsections, we describe the simulations addressing these questions and we summarize our findings. We find that the most contacts occur when the mosquitoes use a crosswind strategy and are allowed sufficient time to either find a host or leave the domain permanently. We conclude that the larger of two unequally-sized groups has a smaller per capita number of contacts. Finally, we find that the number of contacts increases as hosts occupy a larger subregion within the computational domain.

Assessing plume-finding behavior
Intuitively, a crosswind flight strategy should result in a larger number of contacts than upwind or downwind behavior if the plume is straight. This is because mosquitoes close to the plume are more likely to intercept it if their motion is primarily crosswind. But for a meandering plume, it is not obvious if a crosswind flight strategy is more effective than an upwind or downwind strategy. We simulated mosquito behavior in both straight and meandering plumes and recorded the proportion of mosquitoes that found a host and the average time that it took a mosquito to locate a host. We estimated the effectiveness of each plume-finding behavior using these results and found that a crosswind strategy is superior to upwind and downwind strategies, but is less efficient.
The odor plumes were produced by a regular arrangement of nine hosts with a density of 1 host per 1 ft 2 (0.09 m 2 ) located in a single small patch in the center of the square simulation domain. The time series of random velocity fields superposed over the large-scale flow was the same for all simulations. See Table 1 for parameter choices and Eq. (2) for the formula for the meandering plume. An example of the general form of the meandering plume is shown in Fig. 5, alongside a straight plume for comparison. The meandering plume covers more area and achieves a greater width than the straight plume. The outermost contour in both plots is C 0 , the CO 2 sensing threshold of the mosquito. The other contours are equally divided between C 0 and the maximum concentration within each plume. The meandering plume has a  Table 1 for the output variables P (proportion of mosquitoes that found a host) and T avg (the average time to a contact). The second column is the output variable with its baseline value at the parameter set given in Table 1. The third column lists the input parameters associated with the highest sensitivity indices. The fourth and fifth columns are the sensitivity index SI and its error from Eqs. (6)  maximum concentration that is approximately three times higher than that of the straight plume. Since the CO 2 sources are the same in both simulations, the concentration difference is due solely to the differing velocity fields. The results in Table 3 show the proportion of mosquitoes that found a host, P, and the average time to locate a host, T avg . The column labeled P 35 gives the value of P after only 35 s of host seeking. We make the following observations from the table.
N Given unlimited time, the proportion of mosquitoes that found a host was larger in the meandering plume for all plumefinding behaviors, but it took much longer on average to reach a host.
N In both plumes, a larger percentage of mosquitoes found a host using the crosswind strategy, but they took substantially longer to do so.
N If the mosquitoes were time-limited to 35 s, then results for upwind and downwind behaviors were essentially unchanged, but crosswind behavior went from being the most effective strategy to being the least effective in the meandering plume.

Assessing host distribution
We explored the number of contacts that occur in two host groups of equal density but unequal number. An equal per capita contact rate across both host groups would indicate that distributing the hosts into two separate patches has no effect. However we found unequal per capita rates between groups, indicating that a mosquito is less likely to detect a large group of hosts than it is to detect the same number of hosts distributed in smaller groups or as individuals. We also found unequal numbers of contacts between two unequally-sized groups, indicating that one group was easier to find than the other.
We performed a set of simulations in which we held the total number of hosts constant in the domain (10 birds), but split them between two groups in pairs of 9 and 1, 8 and 2, etc., down to 5 and 5 hosts per group. Fig. 2 shows an example of the 7-3 distribution. All other parameters were the same as in the previous section for the straight plume. The two straight plumes from the host groups were well separated from each other and from the left and right domain edges. The exact positions of the hosts were allowed to vary from one simulation to the next, and this affected the plume shape and the internal distribution of CO 2 over the plume.
To compare results between groups we plotted the ratio of the mosquito-host contacts in the smaller group divided by the mosquito-host contacts in the larger group (S/L) against the ratio of the number of hosts in the smaller group over the larger group (see Fig. 6). Each point in the error bar plot shows the mean and standard deviation over 150 simulations for each S/L ratio. The diagonal line in Fig. 6 denotes the case when the per capita contact rates are the same between groups. A value of 1 on the y {axis represents the case when the number of contacts was the same between both groups. When there are 5 hosts in both groups the number of contacts is the same, confirming that there is no leftright bias in the velocity field. When the group sizes are unequal, the per capita contact rates are higher in the small group, but the total number of contacts is higher in the large group (except for the nearly equal 6-4 distribution). This occurs in the region between the lines y~1 and y~x. The crosswind strategy is closest to having equal numbers of contacts in both host groups, which corresponds to equal ease in locating either group.

Assessing changing host density
In this section, we consider the effect of varying host density given a constant number of hosts. Intuitively, the size of the region where the hosts are congregated will affect the number of mosquito-host contacts for two reasons. First, a larger host area is more likely to be found by mosquitoes; second, the spatial arrangement of the hosts affects the size and shape of the odor plume they generate. We performed a set of simulations in which 10 hosts were distributed in a single patch in the center of the domain. Patch area was varied from 10-80 ft 2 (1-7.4 m 2 ) with a corresponding host density varying from 1-8 ft 2 per host (0.1-0.74 m 2 ). Our CO 2 simulation region was 1076 ft 2 (100 m 2 ), so that the patch occupied less than 10% of the simulation region. For each host density, we ran 150-450 simulations to average the effects of host position and individual mosquito choices, with more simulations performed for larger patch areas. The third and fourth columns, P and std dev, are the average and standard deviation of the proportion of mosquitoes finding a host taken over 15 simulations of the same plumes with stochastic mosquito behavior. Each simulation is sufficiently long to ensure that all the mosquitoes either find a host or leave the domain. T avg is the average time to a host taken over all simulations. The associated standard deviation taken over the means of the simulations. The final column recalculates P assuming that the simulation halts after 35 s of host-seeking. doi:10.1371/journal.pcbi.1002500.t003 The proportion of mosquitoes P that made contact in the group is shown by the solid markers in Fig. 7, with the bars corresponding to + one standard deviation. All three plumefinding behaviors exhibit a slight upward trend in P with increasing host patch area. Crosswind plume finding resulted in the most contacts, while upwind and downwind plume finding exhibit nearly identical results. In Fig. 7, P is plotted against the percentage b|100%, where b is the side length of the square host patch divided by the square computational domain side length. b is a dimensionless ratio that relates a length scale important to the host distribution to one that is important to the mosquito distribution. As b increases, the hosts are spread over a larger area and the host density decreases.
We approximated the results by the line P&P 0 zb(1{P 0 ), where P 0 is the proportion of contacts made in the hypothetical case where the area of the host patch is zero. P 0 depends on plume-finding behavior, wind velocity, number of hosts, and other simulation parameters. We performed a least squares fit to find P 0 for each plume-finding behavior (P 0,uw &0:16, P 0,dw &0:17, and P 0,cw &0:33) and the lines are shown in Fig. 7. In the Discussion, we present a first attempt to link this linear approximation from our agent-based model to a contact rate that has been used in standard (nonspatial) epidemiology models.

Discussion
We developed an agent-based/continuum model to explore the effect of behavioral decisions and spatial heterogeneity on the contact rate between mosquito vectors and bird hosts and used it to address three issues of potential interest to researchers in epidemiology and vector control.

Assessing plume-finding behavior
Our results generally show that crosswind plume finding most reliably led mosquitoes to a blood meal source. The effectiveness of the crosswind strategy over periods of minutes is compatible with the conclusions of the mark-release-recapture studies of Anopheles gambiae Giles in [34], where the dispersion of recaptured mosquitoes was related primarily to the distribution of human settlements (over time scales of days). We find that this effectiveness was accompanied by a high cost in time, as also seen in [12]. If the success of host location was restricted to occur within 35 s of mosquito release, then crosswind flight was the superior strategy only in a straight odor plume. In a meandering plume under the time limit, up-and downwind searching were better than crosswind flight. These results are similar to the conclusions drawn by Sabelis and Schippers [20], who used a geometric argument to show that crosswind plume finding is less effective when wind direction is highly variable.

Assessing host distribution
When hosts were arranged in two groups of different size, the larger group consistently attracted more mosquitoes regardless of the plume-finding strategy, although the difference was less pronounced for crosswind flight (Fig. 6). At the same time, the smaller group experienced a higher per capita contact rate. These results are consistent with Foppa et al. [7], which reports on indoor experiments of Cx. quinquefasciatus feeding on roosting chickens. They found that the per capita feeding rate on a single chicken was about 4.27 times higher than that on a chicken in a group of nine. We find similar mean per capita ratios (4.3 and 4.4) for the upwind and downwind plume-finding behaviors when ten hosts are split into a group of nine and a singleton. The close match is surprising since our simulations included wind, whereas the experiments were conducted indoors. However, the uncertainty in the simulations and experiments is high.
The results indicate that it is advantageous for birds to roost together in larger groups on this spatial scale, because on average they will receive fewer bites. This phenomenon of attack abatement is well-known in the literature on predator-prey interactions [35], and likely has two causes in the situation considered here. First, the odor plume exuded by a group does not grow linearly with the number of hosts, because the hosts are clustered together. This leads to fewer mosquitoes locating the group than would find the same number of well-spaced individuals (an avoidance effect). Secondly, mosquitoes only need a fixed amount of blood, and so they will not attack additional individuals even if they are available. This is a dilution effect, also seen in the reduction of groups at risk for malaria resulting from urbanization [36]. The resulting contact heterogeneity arising from attack abatement could have important implications for transmission dynamics.

Assessing changing host density
Standard models of mosquito-borne transmission assume that the mosquito contact rate on one host is inversely proportional to the numbers of hosts (reviewed, e.g. by [37]). This is likely true when hosts are so abundant that a mosquito will always be able to locate one. Our simulations indicate that the probability that a mosquito will locate a host is largely determined by the shape of the odor plume (see Table 3). However, the shape of an odor plume is difficult to predict, since it depends on local wind velocity and turbulence due to landscape features. We therefore propose an approximation that does not depend on plume shape, and is instead based on b, which can be viewed as a patchiness parameter (Fig. 7). This allows us to model contact rates when hosts are variably distributed over patches within a larger domain more realistically even if exact features of the relevant odor plumes are unknown.
We derive a patchiness-driven contact rate model by starting with the contact rate for malaria transmission from Chitnis et al. Figure 7. The proportion of mosquitoes finding a host as a function of host density. The simulation results are given by the solid markers with + one standard deviation (UW = upwind, DW = downwind, CW = crosswind). The lines are linear fits to the data using the formula P~P 0 zb(1{P 0 ), where P 0 is a free parameter and b is the length of a side of the square host patch divided by the length of a side of the simulation region. As b increases, the hosts spread out (i.e., decrease in density) and the number of mosquitoes locating a host increases. The x{axis is b labeled as a percentage. doi:10.1371/journal.pcbi.1002500.g007 [38] which applies specifically in the case of homogeneous mixing in a fixed area without considering host-seeking mechanisms. They propose a vector-host contact rate of where s v is the number of contacts each mosquito wants per unit time, s h is the maximum number of contacts a host can receive per unit time, N v is the total number of vectors and N h is the total number of hosts. This implies that s v N v is the approximate contact rate when the mosquito population is small since every mosquito in a small population is expected to locate a host when the populations are homogeneously mixed. We seek to modify this expression to account for host patch area, as depicted in Fig. 1, and mosquito behavior. When the hosts are limited to a small patch of area within a larger region, we make two modifications to the formula above. First, not every mosquito in a small population can be expected to find a host, since not all of the mosquitoes will come into contact with the odor plume. For this reason, we replace the contact rate with bs v N v , where b represents the ratio of the diameter of the patch where the hosts are located to the diameter of the CO 2 simulation region. The second modification comes from the fact that even when the patch size becomes small (and b&0), the number of contacts must be nonzero since even a tiny area produces a plume that mosquitoes can follow to a host. We want our modified contact rate function to be consistent with the original function in [38] when b~1, which is the case where the host patch is equal to the entire CO 2 simulation region. Based on these observations we propose the modified contact function where P 0 represents the proportion of mosquitoes that find a host as the patch of area containing the hosts becomes infinitesimally small. We note that when the host patches are small (i.e. small b), Eq. (8) can be linearized to read C m &s v N v ½P 0 zb(1{P 0 ), where the factor in brackets, P 0 zb(1{P 0 ), is interpreted as the proportion of mosquitoes that find a host for a given patch-to-region length ratio b. This proportion can be computed from the simulations, as shown by the lines in Fig. 7. Therefore, reasonable assumptions about the patchy distribution of hosts in the domain allow us to devise a quantitatively reasonable estimate of the host-finding probability of mosquitoes for the space and time scales considered here.

Limitations
The length and time scales in the simulations presented here were motivated by the indoor experiments in [7] and a desire to quantify mosquito success near the end of a host-seeking flight. Our conclusions must be validated for spatial domains or time periods that are significantly larger. Furthermore, there are other probability distributions that can be used for mosquito navigational choices (see e.g. [12]) and a variety of initial mosquito spatial arrangements that may affect our conclusions.
Such parameter choices are within the capabilities of the model, and to demonstrate we briefly present a simulation of a growing odor plume in a large domain. Mosquitoes engaging in all three types of plume finding (2000 mosquitoes each) were initially uniformly distributed over a disk of radius 100 m centered on four clusters of ten roosting birds each. The wind meandered according to an incompressible flow appropriate for the large domain size and the crosswind mosquitoes flew in the same direction for a number of decisions chosen from a Lévy distribution having a power law in the tail of p(')*' {2 [39]. The Lévy distribution was centered at 0.7 s, the mean value of T cwd in Table 1, and only positive values were sampled. The mosquitoes were placed in the domain after 150 seconds when the odor plume was approximately 35 meters long. Fig. 8 shows the state of the odor plume after 250 s and 500 s, along with the positions of all nearby mosquitoes. The simulation was run for 1500 s, at which point the odor plume was about 180 m in length and 4.2%, 6.4%, and 26.6% of the upwind, downwind, and crosswind mosquitoes had located a host respectively. The success rate was lower compared to the smaller domain, particularly for the upwind and downwind mosquitoes that rapidly left the plume behind (see Fig. 8, right). The crosswind mosquitoes moved downwind at the same rate the plume did, and therefore had many more opportunities to intercept it. There were still a substantial number of crosswind mosquitoes interacting with the plume when the simulation ceased.
We hope to expand our careful analysis of the smaller domain to a larger domain in the future. In larger domains, it would be interesting to include the effect of mosquito breeding sites and how their locations affect the biting rate of hosts near them in comparison to hosts living farther away. Significant differences in the biting rates in spatial relation to breeding sites have been reported for malaria in [40].

Future directions
There are other potential model extensions that are equally interesting. Disease transmission via mosquito bite, host movement, infection, and demographic processes in both vertebrate hosts and mosquitoes and the dependence of these processes on biotic and abiotic factors could be integrated with the existing model for explicit small-scale modeling of disease spread. This modeling framework is capable of accommodating many further levels of complexity, such as gusting wind, moving hosts, multiple host types, odor-baited traps, variable breathing rate, compound odors, repellents, etc. The challenge will be to identify the components that most strongly affect the behavior of the model system and the underlying reality on which it is based. For example, in the simulations presented here the spread of the host odor is dominated by turbulent convection and diffusion plays a very minor role. It is therefore reasonable to hypothesize that the particular odor cue used by the mosquitoes will have little effect, and that substituting lactic acid for CO 2 (for example) will not impact the resultant contact rates. See [13] for a discussion of turbulent mixing versus molecular diffusion in the context of chemotaxis.
Additional factors could be included in the model in order to capture subtle differences between mosquito species, lighting effects, and other elements not included in the current work. Hosts tend to attract mosquitoes in unequal ways [41]. Differential attractiveness, the emission of different levels of CO 2 by different hosts as well as multiple odor cues can also be introduced into the model to study situations like those presented in [42], where mosquitoes of many species finding humans distributed in huts, or in [43], where mosquitoes were collected inside village houses in Tanzania. These studies report an approximate direct relationship between the number of inhabitants per house and the number of mosquitoes collected. Our model could be used to study which factors influence more strongly this relationship.
From the point of view of controlling the vector population, the model presented here may offer some insights into how the spatial distribution of mosquito traps may affect the overall control. This can be accomplished by replacing the hosts in the model with odor-baited mosquito traps and adjusting appropriate parameters. This was addressed in a non-spatial model by Okumu et al. [10], where homogeneous mixing of hosts was assumed. They discuss the importance of space in the vector-host contact process and indicate that the rate at which an individual host is discovered by an individual vector depends on the distance between hosts and vectors as well as on the size of the odor plumes generated by the hosts. Further, spatial characteristics such as the topography and wind direction are known to be influential in the rates at which individual hosts are found. Our model can explicitly include spatial features to compare strategies of where to place mosquito traps relative to blood-source hosts. According to the studies in [44] the distances at which various species of mosquitoes responded to CO 2 baits by initiating orientation toward the them was 30 meters or less. Therefore, the small model length scales presented here are appropriate for such simulations.