Modeling the average population of La Crosse vectors in Knox County, Tennessee

ABSTRACT La Crosse Virus (LACV) is an arbovirus found in Eastern Appalachia and can cause pediatric encephalitis in prepubescent children. To assess the risk and transmission of this disease, it is particularly important to understand the average population of Aedes mosquitoes, which are the vectors of this virus. We use a deterministic compartmental model to study the effects of environmental factors on the population dynamics of Aedes mosquitoes in the Knox County area. We use locally-collected mosquito population data to adjust our model outputs and find that model transitions are heavily dependent on the fluctuations of both temperature and accumulated precipitation. These findings should be considered for mosquito management in Southern Appalachia, as well as in other regions with slight modifications to our model.


Introduction
is the most frequently diagnosed pediatric arbovirus in the continental United States (Thompson, Kalfayan, & Anslow, 1965). Children under the age of 15 (typically 5-9 years old) infected with LACV may develop symptoms such as headaches, fever, disorientation, and seizures and are then diagnosed with La Crosse encephalitis (LACE) (Erwin et al., 2002). LACE is a significant public health concern in the Midwestern (Calisher, 1994) and Appalachian (Erwin et al., 2002;) regions of the United States, and infection rates often spike in summer and in early fall (Szumlas et al., 1996). The time between when a susceptible human is bitten and when symptoms emerge, or the intrinsic incubation period, ranges from 5 to 15 days (Centers for Disease Control and Prevention 2016). In one study, 12% of hospitalized children with LACE had neurological deficits by time of discharge, and increased levels of behavioural problems were recorded 1-1.5 years after infection; thus, children infected with LACV face the lifelong risk of adverse health outcomes (Erwin et al., 2002).
In the east Tennessee region, where LACV is endemic and is transmitted via the bite of infected Aedes mosquitoes, LACV is maintained in the environment in a zoonotic cycle, between the primary vector Ae. triseriatus, and reservoir hosts such as common gray squirrels (Sciurus carolinensis) and chipmunks (Tamias striatus) (Beaty & Calisher, 1991;Grimstad, 1988). Hypothesized accessory vectors include Ae. albopictus and Ae. japonicus (Gerhardt et al., 2001;Westby, 2015). Although hematophagy occurs, infection with LACV does not have negative health impacts on reservoir hosts (Borucki, Kempf, Blitvich, Blair, & Beaty, 2002). LACV is also maintained within the mosquito population via venereal transmission from male mosquitoes to female mosquitoes (Beaty & Thompson, 1975) and transovarial transmission from an infected female passing LACV to her eggs (Hughes, Gonzalez, Reagan, Blair, & Beaty, 2006). No Known negative health impacts to the offspring (Patrican & DeFoliart, 1985). Eggs infected with LACV retain the virus through their overwintering stage and will ultimately emerge the following season as infected adults (Beaty & Calisher, 1991;Grimstad, 1988).
Thus, understanding the dynamics of the vectors is important for managing LACV and preventing LACE. A mathematical rigorous approach to modelling a population of a species with the life-history divided in age classes has been previously outlined (Gurney, Nisbet, & Lawton, 1983). This work points out that some rates may have a type of delayed effect (Beck-Johnson et al., 2013). Environmental factors, such as temperature and precipitation, are known to influence mosquito population dynamics, and previous studies have shown a relationship between such environmental variables and arboviral infection rates for other vectors and vector-borne diseases (Paull et al., 2017). The different impacts of environmental factors on individual species is not well known and the fine details of how they influence population dynamics have not been well established. Mosquito development from one life stage to another often depends on either temperature, precipitation, or both; environmental factors over an extended period likely play a significant role on lifespan, parity, and biting frequency. The combined effects of rainfall and temperature had an impact on West Nile virus transmission over a single season (Shand et al., 2016); though the specific impacts of temperature and precipitation are likely to vary for LACV compared to West Nile virus because they are transmitted by different mosquito genera. Other environmentally-driven models have been created for Ae. albopictus in a Mediterranean climate (Cailly et al., 2012;Ezanno et al., 2015;Tran et al., 2013).
While LACV is endemic in southern Appalachia, the region does not have an Aedes surveillance programme. In eastern Tennessee where LACV is frequently diagnosed in 20 different counties which vary each year, only one county health department has a mosquito surveillance programme and that programme targets West Nile virus. Most mosquito collections in the region are reactive by the state health department and are in collaboration with local universities (Erwin et al., 2002;Gerhardt et al., 2001;Lambert et al., 2015;Trout Fryxell et al., 2015). If the region had established surveillance, like others do for West Nile virus, then multi-season models of mosquitoes could be developed for the region. To our knowledge, few models exist to study the dynamics of LACV vectors (Bewick, Agusto, Calabrese, Muturi, & Fagan, 2016;Nance, Trout Fryxell, & Lenhart, 2018), furthermore, the potential to use environmental factors as predictors of LACV risk has not been fully realized. As a result, the subsequent benefit of anticipating and responding to changes in LACV risk has not been harnessed. Accordingly, the goal of this project is to use mathematical modelling to develop a model that incorporates temperature and precipitation to explain the dynamics of LACV vectors over a single season in Knox County, TN. This model will use data from a single season from Urquhart et al. study in Knox County in 2013(Urquhart, Paulsen, Moncayo, & Trout Fryxell, 2016.This modelling effort is warranted by the number of LACE cases in the region and the need to establish an Aedes/ LACV surveillance programme; consequently, we begin this modelling effort by modelling the population over a single season. Subsequent developed models will ultimately be used by health departments to prepare for LACV risk periods, and may be used as the basis of similar models for use in other regions experiencing LACV risk. In the next section, we present our model and some details about the corresponding data and its collection. In the third section, we discuss the choice of specific functional forms for rates depending on temperature and precipitation and present those rate functions with parameters estimated from the data from Urquhart et al. (2016). Our numerical results with the estimated rates are shown and are compared with the data. The final section gives a discussion of our work.

Model formulation
We created a model with a system of ordinary differential equations (ODE) to illustrate Aedes abundance over a single season in east Tennessee, specifically for an average Knox County mosquito collection site using field-collected Ae. triseriatus, Ae. albopictus, and Ae. japonicus. All three of these mosquitoes were combined in the model because they are confirmed or hypothesized vectors of LACV, have similar life histories, and were not identifiable at the egg stage; thus, modelling the three species together was justified. Biological and environmental data were used as inputs in the mathematical model, with a particular focus on the impacts of temperature and precipitation in determining mosquito population dynamics. Our model incorporated parameters from local field-collected data obtained from the Trout Fryxell laboratory. This data set was also used in Urquhart's study in Knox County in 2013. Briefly, mosquito traps targeting host-seeking mosquitoes, gravid mosquitoes, and mosquito eggs were placed at eight sites throughout Knox County, TN. Five sites were located at households from 2011-12 positive LACE cases, and the remaining three sites were set within five kilometres of the previous sites, thus selecting for areas likely to have LACV-positive mosquito populations (Urquhart et al., 2016).
The original study by Urquhart et al. (2016) was designed as a trapping efficacy study; consequently, many traps were used at the different sites to identify the best trapping method for each LACV vector. These traps were originally designed to assist collectors in determining the relative mosquito population size. Because they conducted a time-series monitoring effort with a variety of traps targeting different life stages and environmental conditions, we used these data to generate our initial model. We hope to continue these efforts and be able to forecast LACV transmission in the region. The data used here will resemble future surveillance efforts as most surveillance activities are limited to monitoring relative vector abundance; these data are often what is used to drive management decisions.
Daily temperature and precipitation data were obtained from the National Oceanic and Atmospheric Administration's (NOAA) National Weather Service and National Centers for Environmental Information, respectively. Since the amount of standing water present at a collection site greatly influences the oviposition process (Wang, Tang, & Cheke, 2016), we utilize accumulated precipitation prior to mosquito collection instead of daily precipitation in our model. Temperature and accumulated precipitation were recorded over the study season. The daily mean temperature in 2013 from June 29 (calendar day 180) to October 27 (calendar day 300) is plotted in Figure 1(a). The temperature fluctuates until the hottest period of summer, around calendar day 240, after which the temperature trend starts to decrease. Studying over a large range of temperature provides us the opportunity to effectively map the effects of temperature on mosquito dynamics. Precipitation was summed seven days prior to mosquito collection, termed accumulated precipitation, since simulations were most accurate with seven-day accumulation (described in detail below). Accumulated precipitation is highest in the beginning of the study period, varying throughout the following weeks and is displayed in Figure 1 The model diagram in Figure 2 shows the model compartments. We focus our model on female mosquitoes only. Compartments represent mosquito life stages: (E) egg, (I) immature, (H) host-seeking, and (G) gravid. Each compartment represents the average number of mosquitoes per trap area in the corresponding life stage. The egg compartment consists of mosquito life stages from oviposition and embryonating to hatching. The immature compartment consists of the several larval life stages and pupae. The host-seeking compartment consists of adult mosquitoes searching for blood meals. Finally, the gravid compartment consists of mosquitoes ready to oviposit. In the model system below, T(t) represents the temperature at time t and P(t) represents the accumulated precipitation at time t. The differential equations above describe the movement of mosquitoes across life stages. The compartments have fixed mortality rates, with μ E , μ I , μ H , and μ G relative to the concurrent life stage of mosquitoes. Additionally, since the typical mosquito lifespan is spent at a single collection site representing ∼300 metre in radius, no migratory terms were included in our model (Medeiros, Boothe, Roark, & Hamer, 2017). The ODE for the egg compartment contains the term b(P)G, which represents that eggs laid by gravid mosquitoes enter the egg compartment at the oviposition rate of b(P). While male mosquitoes can carry LACV if infected transovarially, only female mosquitoes can infect humans by biting. To represent this, the model does not consider male mosquitoes. Hence, the oviposition rate b(P) was adjusted to represent a one-to-one sex ratio. The transition rate at which eggs hatch and reach the larval stage is represented by the function f E (T, P). The term f E (T, P)E represents the movement from the egg compartment to the immature compartment. Finally, μ E describes the death rate for eggs, and the term μ E E describes the deaths in the egg compartment. The specific structure of the functions T(t) and P(t) will be discussed in the next section.
The immature compartment has entry from the egg compartment at the rate of f E , represented by the term f E (T, P)E. The transition rate at which larvae emerge and mature into host-seeking adults is represented by the function f I ; the term −f I (T, P)I represents the movement of mosquitoes from the immature to the host-seeking compartment. Similar to the ODE for the egg compartment, the corresponding death rate for the immature compartment is denoted by μ I .
Mosquitoes enter the host-seeking compartment at the rate −f I (T, P) from the immature compartment, represented by the term −f I (T, P)I. Once the host-seeking mosquitoes find a blood meal, they move to the gravid stage ready to oviposit. This movement occurs at the rate of d, is represented by the term −dH. Once the gravid mosquitoes have oviposited, they move back to the host-seeking compartment at the rate of γ and restart the process of finding a blood meal. This process is represented by the term γ G. This compartment also has a corresponding death rate of μ H . The gravid compartment has entry from the hostseeking compartment at the rate of d, represented by the term dH. After oviposition, the gravid mosquitoes return to the host-seeking compartment at the rate of γ , represented by the term −γ H. This compartment also has a death rate of μ G .

Results
As evident in Table 1, most of our parameters were estimated from the field-collected data. While b(P), f E , and f I are functions of T and P, the remaining parameters are constants in the model.  (Brady et al., 2013) To include the influence of environmental factors on the mosquito population dynamics, we set the oviposition rate b, the transition rate from eggs to immature mosquitoes f E , and the transition rate from immature mosquitoes to host-seeking adults f I as functions of either temperature or accumulated precipitation. We determine the functional forms for b, f E , and f I based on how mosquito populations change over time. We also assume that temperature and accumulated precipitation affect these three rates differently. Note that there is no evidence that the natural death rates of these populations depend on temperature and precipitation. The transiton rates between the gravid and host-seeking populations have been estimated in Trout Fryxell's lab.
We use ODE45, a higher order Runge-Kutta method in MATLAB (MATLAB, 2017), to numerically solve our system of differential equations. To find functional forms for b, F E , and F I , we tested combinations of polynomials, logistic functions, and Holling's type II functional responses. We also decided how many days to include in accumulated precipitation; we checked accumulation over 1 day, 2 days, and continuing up to 10 days. In this fitting procedure, we searched over parameters for each functional form and over each length of the accumulation. We set the objective functional as the relative error of the MATLAB simulated mosquito populations from the field-collected mosquito populations. Then, to estimate parameters and various functional forms, we minimized this relative error through the Ordinary Least Squares method, where the data were the totals of eggs and the combination of gravid and host-seeking adults on observed days. Due to the small number of gravid mosquitoes (n = 300 total) compared to hostseeking mosquitoes (n = 2672 total), we fit the model to host-seeking mosquitoes. Also, because gravid mosquitoes are only gravid for about 2 days before they become hostseeking mosquitoes and because host-seeking mosquitoes transmit LACV we considered host-seeking mosquitoes as the more important life stage.
We used the MATLAB function fminsearch (a simplex algorithm used for constrained optimization) to narrow the parameter ranges and then used fmincon (an interior point algorithm used for constrained optimization) to estimate the best parameters for each functional form and for accumulation level. The functional forms with their corresponding fitted parameters and length of accumulated precipitation with the smallest relative error were used in simulations. To give some more background on our choices, we discuss below more details of the mechanisms drivng the choices of functional forms.
Since we know that cumulative precipitation four weeks prior to collection may indicate levels of Ae. albopictus populations (Haddow, Gerhardt, Jones, & Odoi, 2009), we hypothesize that with increasing accumulated precipitation, the oviposition rate b(P) continues to increase until the accumulated precipitation gets close to a certain level. After that point, the oviposition rate stays approximately constant. Hence, we used functional forms which act accordingly to our hypothesis; a logistic type function (plus a constant) was the rate chosen for b(P).
Since the transition rates from eggs to immatures and immatures to host-seeking are affected by both precipitation and temperature, (Nance et al., 2018;Wang et al., 2016;Xu et al., 2016), f E and f I are functions of temperature and accumulated precipitation. In our choices, we considered functional forms that had T + P as variables and also those had two separate functions of T and P added separately. We hypothesize that with increasing temperature, the transition rate from eggs to immatures increases and eventually plateaus. For f E , the sum of a Holling's Type II functional responses in E and a logistic function in P (plus a constant) were chosen.
We hypothesize that with increasing temperature and precipitation, the transition rate from immatures to host-seeking increases until a certain point, after which this transition rate decreases with increasing temperature and precipitation (Delatte et al., 2009). A downward facing parabola in the sum T + P was chosen. The selected rate functions for b, f E , and f I are listed below. Initial values were obtained from field-collected data. Due to the under-sampling of the first three weeks, we begin simulations on week 4. At that time, the Trout Fryxell laboratory collected 291 eggs and 31 mature mosquitoes from the field. In this context, both host-seeking and gravid mosquitoes are defined as mature adults. To determine the initial values for the remaining compartments, we assume that the immature population was roughly half as large as the egg population. Additionally, there were slightly more host-seeking mosquitoes than gravid mosquitoes in our data. Thus, simulations were run with E = 291, I = 145, H = 17, and G = 14 as initial conditions. With our estimated parameters and functional forms, we run our simulations.
Our simulations for each one of the compartments in our model are shown in Figure 3. We plot the summer months against the average abundance of mosquitoes per site. The abundance is high throughout the summer; however, the mosquito abundance plummets towards the end of the summer months as expected due to the decrease in temperatures and in accumulated precipitation. In Figure 4, we compare our simulations for adults and eggs against the field-collected data for eggs and adults. To check the accuracy of our simulations, we compute the relative error. We find that the relative error for the simulation for adult compartments is 20%, and the corresponding relative error for egg simulations is 25%.

Discussion
Using data from 2013 collected by the Trout Fryxell laboratory and presented in Urquhart et al. (2016) we chose our functional forms and parameters and length of accumulation to represent mosquito populations at an average sampling location in Knox County, TN. This model demonstrates the value of integrating field-collected data with mathematical modelling, as we have created a model that is both sensitive to environmental conditions and adjusted for specificity to the geographic region in which it will be applied. Furthermore, by modelling a single site/location, our model is more applicable to the field, as a single site more closely matches the extent to which humans are exposed to mosquitoes and their associated viruses. In fact, our model can be easily adopted in different regions across the world.
In application, the results of the model indicate a few key conclusions relative to pest management and public health efforts. The impact of accumulated precipitation on the model reflects the importance of drainage and the need to prevent the formation of standing water in areas frequented by those likely to contract LACV. Our study provides evidence that water from precipitation plays an important role in the life cycle of mosquitoes. Additionally, the results of this study indicate that modelling can be a useful part of integrated pest management, as informative campaigns and strategic pesticide use before and during predicted peak times can generate public awareness, minimize undue concern, and allow for efficient allocation of pest management resources.
A number of errors reduce the accuracy of the model's fit, though these errors do not disqualify the validity of the model overall. The use of a single season's data creates a liability for overfitting the model to conditions unique to that single season. Weather variations between the same seasons of different years could indicate the need to do future work incorporating multi-year seasonal data in a model. Variations between generalized weather data and microclimates at sampling locations may also account for error in the fit of the model (Murdock, Evans, McClanahan, Miagowicz, & Tesla, 2017). Additionally, the absence of environmental terms aside from temperature and precipitation may create limitations, when the impact of other weather events are substantial. Using temperature and precipitation in the model allows incorporation of the environmental factors most influential on mosquito population dynamics. Previously, Paull et al. (2017) identified drivers for West Nile virus epidemics in North America. In that paper, drought was a driver for WNV and in Tennessee it was a combination of increasing temperatures and decreased precipitation that was associated with West Nile virus-infected mosquitoes (Paull et al., 2017). However, other environmental factors, such as humidity, wind-speed, local flora, and amount of shade also likely play a role in mosquito population and activity levels. For this reason, incorporating terms for other weather factors may be a useful future addition. Finally, the error in fitting egg abundance was much greater than that of the adult mosquito population. Adult mosquitoes are responsible for transmission of LACV, and we therefore prioritize creating a model that accurately fits the population dynamics of adult mosquitoes in general; however, future improvements in the model's egg abundance calculations may result in more accurate predictions of egg abundance and adult mosquito population dynamics.

Acknowledgments
We want to thank Cassandra Urquhart and Dave Paulsen with University of Tennessee Institute of Agriculture's Department of Entomology and Plant Pathology for collecting and identifying the mosquitoes. We also thank Brian Hardison for useful comments.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was conducted, in part, during the 2017 Summer Research Experience (SRE) for undergraduates programme at the National Institute for Mathematical and Biological Synthesis, an Institute sponsored by the National Science Foundation through NSF Award #DBI-1300426, with additional support from The University of Tennessee, Knoxville. Entomology work was completed using funding from University of Tennessee Hatch Project TEN00433 and Abelardo Moncayo with the Tennessee Department of Health.

Data Access
The data that support the findings of this study are available from Rebecca Trout Fryxell, RFryx-ell@utk.edu, upon reasonable request.