Application of reaction-diffusion equations for modeling human and breeding site attraction movement behavior of Aedes aegypti mosquito

1.   Introduction

The modeling of dispersal through the use of reaction-diffusion equations has a long history. As early as 1937, it was shown that a diffusion model with logistic growth terms can generate traveling waves ^[1]. From this seminal work sprang numerous mathematical publications on the existence of traveling wave solutions, see, for example, Hadeler and Rothe ^[2] and the review articles of Volpert and Petrovskii ^[3] and Wang ^[4]. Coupled nonlinear reaction-diffusion equations also have a long tradition of application to model the spread of diseases. The classical SIR epidemiological model has been applied in combination with diffusion to generate epizootic waves, as has been modeled for rabies dispersal in England by Murray and Seward ^[5] and Murray et al. ^[6]. There exist many publications on the dispersal of insect-related disease vectors within the frame of reaction-diffusion equations. Maidana and Yang ^[7] analyzed the spread of dengue disease via the disease vector Aedes aegypti and developed criteria for the onset of traveling waves in one dimension. Almeida et al. ^[8] evaluated the use of insecticides and male sterility induction as a means of blocking traveling waves from the mathematical point of view in a one-dimensional setting, providing also a theoretical proof of the efficiency of the sterile insect technique (SIT). Also, in a one-dimensional setting, Anguelov et al. ^[9] proved that bistable traveling wave solutions exist if sterile males are introduced; they illustrated their theorem via numerical investigations. In other papers, the focus is on dispersal in realistic two-dimensional spatial settings. Lutambi et al. ^[10] discretized the space by introducing heterogeneous hexagonal patches. For each patch, a system of ordinary differential equations for the mosquito dynamics was setup, and it allowed for movement between neighboring patches. Yamashita et al. ^[11] carried out two-dimensional numerical simulations of a reaction-diffusion model for Aedes aegypti in a residential environment consisting of streets, blocks and a beach by using a finite-volume method. Knight ^[12] developed a model of the mosquito-mangrove basin ecosystem, detailing the habitat of the saltwater mosquito Aedes vigilax of a mangrove forest in Australia. Richter et al. ^[13] modeled the invasion of a new species by coupling geographic information system data on temperature using a finite-element tool.

Recent experimental studies have shown that the heat emanating from humans ^[14] and olfactory agents from breeding sites ^[15] control mosquito movement patterns. It was found that host seeking by mosquitoes and other pathogen-spreading insects relies on the detection of host-associated cues, including body heat ^{[14,16,17,18,19]}. Heat strongly stimulates mosquito blood feeding, as heat seeking is a part of a multimodal host-seeking program followed by mosquitoes in which body heat serves as an important cue of the host ^{[14,17,20,21]}. After a blood meal from a host, female mosquitoes look for suitable oviposition sites to guarantee success of the offspring. In order to remember its breeding site, mosquitoes need to adapt by learning to respond to odor cues from larval, pupal or early adult vectors ^[22]. Aedes aegypti mosquitoes are capable of information retention, which is acquired either during larval development, hatching or foraging, and that new information could influence their behavior ^[15]. Learning from its own experience could influence the vector behavior, such as the potential preferences in vector-host interactions ^[23,24].

Reaction-diffusion equations with simple diffusion terms do not reflect realistic dispersal patterns of mosquitoes, which are driven by an attraction to prey (humans) and the attraction of breeding sites. The application of reaction-diffusion equations for the dispersal of mosquitoes (and other species) therefore faces two challenges:

1) Coupling the governing equations with geographic information on landscape heterogeneity e.g., types of human settlements and possible breeding sites.

2) Accounting for animal behavior in the equations.

It was the purpose of this study to model and analyze mosquito movement at small scales, including attractant behavior, and to model the dispersal of mosquitoes and spread of diseases in a human population at the landscape scale. To this end, the governing equations were implemented in the finite-element tool COMSOL Multiphysics. For the landscape scale studies, landscape data extracted from Sentinel satellite images were imported using a geographic information system.

The chosen area for this model application is a suburban zone of a metropolitan area in the south of Vietnam. In Vietnam, mosquitoes and mosquito-borne diseases pose severe problems. They occur in most provinces and cities across the country, but are more common in the southern region, especially in accelerating urbanization zones such as the southern area of Ho Chi Minh (HCM) city (Figure 1). Therefore, for our study area, we chose a section of the urbanization zone in the south of HCM city; it has a mixture of vegetation, water bodies, rural settlements and urban residential areas. The study area includes District 7, Nha Be District and the Binh Khanh Ward (Figure 1). The area belongs to the downstream area of the Dong Nai and Saigon river system. This study area was previously covered by a mangrove forest and has slowly been urbanized since 1972. The urbanization has been accelerating since 1990.

The target mosquito species of this study is the Aedes aegypti species, which is involved in the transmission of numerous viruses, including dengue, chikungunya and Zika. There are many studies on the mosquito habitat and life cycle in Vietnam, and the earliest work could be the work of Herbert et al. ^[25], who conducted mosquito surveys and analyzed their distribution. Later, Huber et al. ^[26] studied the ecology of Aedes aegypti, genetic differentiation, variability in competence as a vector for dengue virus and resistance to insecticides; this was to assess the role of the vector in the changing pattern of the disease in Vietnam. Jeffery et al. ^[27] collected and classified water containers from large tanks to small jars in the center of Vietnam to determine if Aedes aegypti populations were spatially and temporally homogeneous.

The objectives of the study were to perform a numerical investigation that analyzed models of the dynamics of mosquito populations that are driven by an attraction to prey (humans) and the attraction of breeding sites, as well as of the spread of vector-borne diseases in urban-wetland areas.

2.   Materials and methods

2.1.   Mathematical models

2.1.1.   Basic model

Notations

State variables

A: mosquitoes (aquatic phase)

M: mosquitoes winged phase (general)

M_s: mosquitoes (not vector of virus)

M_i: mosquitoes (vector of virus)

M₁: mosquitoes (searching for blood meal)

M₂: mosquitoes (searching for breeding site)

M_ij: mosquitoes in search mode i and infectious mode j

H_s: non-infected humans

H_i: infected humans

H_r: removed (immune) humans

Parameters

c(x, y): environmental capacity (aquatic phase)

D: diffusion coefficient

$\phi$: oviposition rate

γ: hatching rate

${\mu }_{A}$: mortality rate aquatic phase

${\mu }_{M}$: mortality rate winged phase

${\mu }_{s}$: mortality rate due to insecticide spraying

d: width of spraying corridor

y₀: position of spraying corridor

r: survival of hatched larvae

ß: contact rate for mosquito-susceptible humans (small-scale model)

ß₁: contact rate for not virus bearing mosquitoes with infectious humans

ß₂: contact rate for infectious mosquito-susceptible humans

${\tau }^{-1}$: infectious period

${\omega }_{12}:$ rate of change between mosquito states

$v$: conservative flux convection coefficient

The general form of a system of reaction-diffusion equations for biological populations with a density y_i is given by

where the spatial operator has the form

with the conservative flux coefficient

which allows for movement in the direction of a gradient of an attractant c. The reaction term f_i accounts for the population dynamic and genetic processes ^[28]. Note that the system may also include components that are stationary, as is the case for a stage-structured population.

For a mosquito population, the basic model consists of the aquatic phase (A) and the winged phase (M) ^[7]. The capacity term varies in space.

The winged phase is coupled to the aquatic phase via the hatching term rγA, as follows:

In Eq (5), the first through fourth terms respectively represent the random dispersal, vectored dispersal, hatching and mortality rate.

2.1.2.   Dispersal patterns at small scales

This model takes into account the blood-feeding status of female mosquitoes, differing with respect to dispersal patterns. Unsaturated mosquitoes hover around their prey, attracted by the body heat. It was recently found that a receptor (IR21a) is activated when mosquitoes move toward cooler temperatures ^[14]. After the blood meal, mosquitoes search for breeding sites for oviposition, guided by olfactory preferences ^[15]. The different movement patterns are expressed by chemotaxis terms; for humans (in terms of temperature), the term is

for breeding sites (in terms of olfactory agents), the term is

The chemotaxis fields are generated by a linear diffusion process with humans and breeding sites as sources and boundary absorption conditions.

Mosquitoes become vectors of disease if they bite infectious humans. Therefore, in the model, the population is differentiated into susceptible (Index s) and infected (Index i) subpopulations. The system is solved over a region G, which includes breeding sites denoted as B, human settlements and vegetation.

Mosquitoes in blood-seeking mode that have hatched at a breeding site ($r\gamma A$) disperse and seek humans (term: ${v}_{T}{M}_{1s}$). After biting either non-infected or infected humans, they switch into the breeding-site seeking mode ($\beta \left({H}_{s}+{H}_{i}\right){M}_{1s}$). The last term represents the switching following oviposition.

Infected mosquitoes originate from encounters with infected humans and do not develop from eggs. Otherwise, their behavior is the same as that of the non-infected type.

Non-infected mosquitoes that have become saturated as a result of biting humans at a rate $\beta {H}_{s}{M}_{1s}$ are attracted to breeding sites (term: ${v}_{B}{M}_{2s}$) and switch back to the blood-seeking mode after oviposition (last term).

Saturated mosquitoes that have become infected as a result of biting infectious humans at a rate $\beta {H}_{i}{M}_{1s}$ are attracted to breeding sites (term: ${v}_{B}{M}_{2i}$) and switch back to the blood-seeking mode after oviposition (last term), where

where $B\subset G$.

2.1.3.   Dispersal patterns at large scales

At the landscape scale (areas of several square kilometers), the movement from local breeding sites to humans is neglected. The objective of models at this scale is to simulate the spread of the disease in consideration of mosquito densities and the densities of human populations. To this end, capacity terms are formulated to reflect the availability of typical breeding sites. Two types of breeding sites are distinguished according to different carrying capacities: water bodies outside human settlements, such as ditches, or small ponds with low carrying capacities and small water bodies in plastic containers, mud dishes, mud pots, ditches and discarded tires with high capacities ^[29]. The latter breeding sites are connected to the types of human settlements. Here, we distinguish rural settlements with lots of such breeding sites and cities with a highly developed infrastructure and not much littering. The governing equations are based on the work of Maidana and Yang ^[7], modifications to the chemotaxis terms and nonlinear dispersal.

Note that parameters such as the environmental capacity or contact rates are dependent on the infrastructure of the study area. The parameter values are shown in Table 1.

2.1.4.   Control management

Release of sterile males.

The effect of this method is the decrease of fertile mating by the factor f = $\frac{{M}_{s}+{M}_{i}}{{M}_{s}+{M}_{i}+S}$, where S denotes the density of sterile males under the condition of a sex ratio of 1:1.

The number of eggs is therefore reduced by this factor and Eq (13) is modified as follows.

The spread of sterile mosquitoes is modeled by a diffusion equation with a mortality term and a source term, which summarize the release rate R_i at each location $\left({x}_{i}, {y}_{i}\right)$.

Spraying in a corridor.

A spraying corridor parallel to the x axis is defined by the y position y₀, the width d of the corridor and the form parameter n, which determines the reduction in pesticide action. The spraying introduces an additional mortality rate.

where ${\mu }_{s}$ denotes the maximum mortality, which decreases with distance from the centre of the corridor y₀.

2.1.5.   Boundary conditions

For the small-scale model, Neumann conditions were applied to the walls of the houses and the zero Dirichlet condition was applied to the outer boundaries. For the large-scale model, the zero Dirichlet condition was applied to the boundaries of the region. The results are only slightly influenced if the outer boundary conditions are Neumann.

2.2.   Methods

The governing equations were implemented using the finite-element tool COMSOL Multiphysics version 6.0 (www.comsol.de). This tool is well adapted for the combination with geographical data. For time-dependent problems, the variable order backward differentiation formula solver in COMSOL was used. For the simulation of dispersal on the landscape scale, COMSOL was linked to a geographic information system that provided population density data and the landscape geometry. Geo-referenced data for the study area were imported from the resultant raster maps for vegetation density level, water body coverage and mean population density distribution.

This information was extracted from Sentinel-2 imagery. Sentinel-2 is a high-spatial-resolution (10, 20, 60 m), high-temporal-resolution (5 days), multispectral (13 bands) imaging satellite carrying a multispectral imager (MSI) ^[30]. The level 1C data of an image acquired on January 17, 2021 were downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu/). These data were pre-processed using a SNAP Desktop platform (Version 8.0.0; SNAP 2020) with the Sentinel-2 Toolbox. This Level 1C image with top-of-atmosphere (TOA) reflectance was processed to a Level-2A bottom-of-atmosphere (BOA) product using the European Space Agency Sentinel Application Platform (ESA-SNAP) Sentinel-2 Toolbox with the additional Sen2Cor plug-in for atmospheric correction. Sentinel-2 bands 5, 6, 7, 8a, 11 and 12 (all are short-wave infrared) were resampled to 10 m × 10 m from their native resolution of 20 m × 20 m to match the resolution of the other bands. These processed bands were used to extract the levels of vegetation cover, the water body and the urban-nonurban areas. The population density was obtained from the Worldpop (www.worldpop.org) global high-resolution population denominators project, which provides gridded population counts at 100-m spatial resolution (3 arc-seconds).

The data were interpolated within COMSOL by using the two-dimensional linear interpolation option in the function menu. The implementation scheme is shown in Figure 2. The environmental capacity of mosquitoes is related to the residential environment, which has been characterized by scores ranging from 0 to 3. In a rural setting (Score 2), the preferred breeding sites of Aedes aegypti, like small water containers, pots or tires, are more common than in dense urban areas (Score 3). Scores 0 and 1 are given to water bodies and vegetation. For numerical reasons a smooth curve was created to avoid abrupt changes in the parameters following transitions of infrastructural features (Figure 3).

3.   Results

We applied this model to investigate the movements of mosquito populations and the transmission of dengue disease to humans by infected vectors. Two scales were setup for the geometry settings. At the fine scale, this model can predict how mosquito populations move between the houses and breeding sites and how the humans are infected by the disease in the area. At the landscape scale, the model can simulate the dispersal of infectious mosquitos and humans and the effects of control measures.

3.1.   Small scale

Here, we have constructed a typical rural setting consisting of small houses and different breeding sites. Between the houses (rectangular shapes), there is a dich, near the houses there are small breeding sites representing containers, pots, and other mosquito-preferred sites. The small circles represent humans and their range of movement. In the model, breeding sites and humans constitute the source of the respective attractants.

The circles in the houses represent humans and their range of movement. In the model, breeding sites and humans constitute the source of the respective attractants. Figure 4 shows the velocity fields of mosquitoes in the non-saturated and saturated states. Due to the chemotaxis effects, one can clearly distinguish between different movement patterns.

Along with mosquito movement, an infection is spread if infectious humans are present. This is shown in the simulation results in Figure 5. The small circles in the houses represent humans and their radius of interaction. At the beginning, one human outside of a house is infected (red circle, Figure 5(a)). Note that humans are fixed in space, i.e., there are no diffusion terms in the governing equations. After 3 days, a portion of the humans are infected (circles with colors, Figure 5(b)); after 5 days, all humans are infected (Figure 5(c)). The movement pattern is determined by the exchange between breeding sites and humans (Figure 5(d)).

3.2.   Large scale

The following simulations show the attractive effects in a real landscape. The dispersal of mosquito populations has been followed from a focal point in the area. Different types of dispersals were observed, including the free dispersal of mosquitoes under the conditions of no treatment regime and the spread of mosquitoes and disease transmission with the aforementioned control management.

3.2.1.   Uncontrolled dispersal

The governing equations, presented as Eqs (13)–(18), are solved for the study area. Here, the detailed movement between breeding and feeding sites is not considered; however, the directed movement to humans is maintained. Figure 6 shows the movement pattern of mosquitoes for a section of the study area.

As was to be expected from the many theoretical studies mentioned in the introduction, traveling waves occur if certain parameter constellations are given. Figure 7 shows a 3D representation of a traveling wave of infected humans, as obtained for a homogeneous domain. In a nonhomogeneous environment, traveling wave fronts are distorted as shown in Figure 8.

Figure 8 shows the wave of the spread of infected humans. Distorted wave fronts reflect the heterogeneity of the landscape. Note that the underlying heterogeneity of the resources and human population density are given in Figure 1.

Accordingly, the density of the susceptible population is decreasing in a complementary manner. Figure 9(a)–(d) show the spatiotemporal evolution of the susceptible population from Day 0 to Day 250. One can clearly see that the disease progresses like a wave with an irregular wavefront, reflecting the landscape heterogeneity.

It is also interesting to look at the time evolution of the epidemic at single points (Figure 10). One can recognize the typical time courses of the classical SIR model, irrespective of the time delay due to the spatial separation. At any given time, individual humans will be in one of three states, namely, S (susceptible), I (infectious) or R (recovered/immune). The infected mosquitoes in the model do not have a recovered/immune state, as it has been assumed that mosquitoes remain infected until they die. One can see that dengue transmission peaks travel along with the traveling wave fronts.

3.2.2.   Control management

An understanding of the movement behavior of mosquitoes and their preferable breeding sites enables the assessment of management strategies (method, frequency, area). Here, we show some simulations of control measures for illustrative purposes. Management methods were incorporated into the model, and they include conventional means of vector control, i.e., insecticide spraying and the use of sterile males to decrease the frequency of fertile mating. One question to answer was as follows: how can a traveling wave be stopped by spraying? In the simulation, a corridor of intensive spraying was implemented. This was achieved by introducing a high mortality rate to a selected strip. The simulations show that the success of this measure depends on the breadth of the spraying zone, as shown in Figure 11.

If the breadth of the corridor is too small, the traveling waves are delayed but are able to persist. So, there exists a threshold value for this parameter, which is related to the minimal viable population size, as discussed by Richter et al. ^[13]. The sterile-male-release technique has recently seen a comeback ^[31]. The effect of this method is a decrease in fertile mating. This control measure has been incorporated into the model in a straightforward way by simply adding a further equation for sterile males and the locations of release (cf. Eqs (19) and (20)). The success of this measure depends both on the release rate and the population density at the location of the release sites. As was proved in the paper written by Anguelov et al. ^[9], the release of sterile males leads to a bistable reaction-diffusion equation. They found that "the dynamics of our system is driven by a sterile insect technique (SIT)-threshold number above which the SIT control becomes effective and drives the system to elimination". For the parameters of our model, this threshold was found to be higher than the environmental capacity. However, due to the property of bistability, elimination can be reached at lower release rates if the insect population is low. This behavior was also immanent in our model equations. Here, we show only one result of many simulations runs.

Figure 12 shows the density of sterile males released at four locations and their effects on the mosquito population.

4.   Discussion

For heterogeneous landscapes, there are quite a few approaches wherein the region is subdivided into patches, each of which with different environmental variables ^[10,32]; the ordinary differential equations are solved for each patch, allowing for migration between patches. Such approaches are similar to reaction-diffusion equations.

Although the theoretical basis of animal dispersal and the dispersal of epidemics in the form of reaction-diffusion equations was established decades ago, such approaches still constitute a powerful tool for realistic simulations, especially when finite-element methods are coupled with geographic information. Furthermore, the underlying mathematical theory has been well elaborated and allows one to decide whether the problems of the solution of the initial value boundary problems are well posed; it also guides us in the understanding and evaluation of the numerically obtained results. Our simulation experiments show that complex movement patterns of mosquitoes can be well incorporated into reaction-diffusion equations. The model takes into account the implication that female mosquitoes retain information from the larval and pupal environment from which they emerge as adult. This was discussed in ^[33,34]. The dissociation between the oviposition site and remote food source (human blood) requires the ability for long- or mid-range orientation, and it would necessitate "some form of learnt behavior" ^[35,36]. Breeding site preferences determine the movement patterns. Gaburro et al. ^[15] found that semiochemicals derived from Aedes aegypti larval rearing affect the future female choice for oviposition. If infected by the dengue virus, preferences for special breeding sites might be lost. Olfactory learning processes were altered during dengue infection. They hypothesized that dengue virus infection alters gene expression in the mosquito's head and is associated with a loss of olfactory preferences, possibly modifying the oviposition site choice of female mosquitoes.

Reaction-diffusion models of the kind presented here, when being coupled with breeding site preferences and driving variables such as the diurnal rhythms of mosquito activity and seasonal changes in moisture and temperature, show that predator–prey interactions can become a powerful tool in predicting infected mosquito behavior and the transmission of mosquito-borne disease at several scales. For long-term predictions, it is also important to consider the evolution of resistance against pesticides.

For future development and application of the model, experiments and observational studies need to be carried out to contribute to model parameters, such as the transition rate between mosquito states, and parameters related to the dispersal and strength of attractants. A promising method for parameter identification is the combination of models and Bayesian statistics, as has been recently applied to the spread of Aedes aegypti in urban areas ^[37]. Once the model is fully parameterized, it can be used for the realistic planning of various control measures, such as the reduction of breeding sites and the reduction of contact rates via application of repellents and biological and genetic control measures.

5.   Conclusions

Although our mathematical model for the movement patterns of mosquitoes, which is based on the mode of selection between blood-meal sites and breeding sites, is just at the theoretical stage, it has potential for real-world applications. When being coupled with other processes and observational data, the model can help to (ⅰ) assess the effects of fast land-use changes and sanitary conditions of an area from the perspective of the dynamics and life cycle of mosquito population, (ⅱ) assess the effects of the heterogeneous characteristics of the landscape and of the mosquito population on the transmission of the vector-borne disease to humans and (ⅲ) assess different landscape designs and the treatment methods to prevent the transmission of mosquito-borne disease. Future work should focus on updating the methods to retrieve plausible parameter sets for the model.

Acknowledgements

The authors wish to thank the anonymous reviewer for his thoroughly reading the manuscript and for many helpful suggestions.

Conflict of interest

The authors declare that there is no conflict of interest.