Creating an Agent-Based Model to Explore the Spread of Zika Virus
The sweeping spread of Zika Virus and its severe symptoms make it a serious threat to public health. However, its records have not become well-kept until recently. The lapses in case reporting have made it impossible to model Zika’s spread by traditional methods, such as differential equation-based models, because they require aggregate data that are not presently available.
We applied an alternative modeling technique and created an agent-based model that is dependent on knowledge of mosquito behaviors and ecological data. Such information is easier to acquire and validate compared to data of the spread patterns of Zika. Our model also enables performing what-if analyses on Zika’s spread under several conditions, predicting and producing results which indicate a causal relationship between outbreak severity and population density. We also show that future work can be focused on refining the behaviors of mosquito agents and enhancing the modularity and scalability of running experimental trials with our model.
Zika Virus (Flavivirus Zika, ZIKV) is an international epidemic. Though it was originally isolated from the blood of a sentinel rhesus macaque residing in Uganda’s Zika Forest in 1947, it has proven its ability to spread in humans lately . As of this writing, the Center for Disease Control (CDC) reports that Zika is being actively transmitted in 62 countries worldwide, amounting to 37,168 cases since the start of 2015 in the U.S. and its territories alone [2,3].
Numerous serious health risks beyond Zika’s placid symptoms—mild fever, rash, joint pain, muscle pain and conjunctivitis—have been noted since the past decade’s outbreaks . One rare yet dangerous side-effect associated with ZIKV after a 2014 outbreak in French Polynesia is Guillain-Barre syndrome, an autoimmune disorder that causes acute and chronic flaccid paralysis [5, 6]. Several studies have also observed a correlation between infant microcephaly and the presence of ZIKV in mothers during the first trimester of fetal development [7, 8]. These high-impact risk factors, which are especially dangerous for young women, combined with Zika’s increasing ubiquity have prompted such agencies as the National Science Foundation to solicit research addressing the ecological transmission dynamics of the virus .
Many mathematical modeling templates have been conventionally used to delineate the spread patterns of epidemic diseases, the Susceptible-Infected-Recovered (SIR) model being among the most commonly adapted of these. In the SIR model, variables such as the number of susceptible, infected, and recovered individuals and the infected, susceptible, and recovered fractions of the total population were joined by the variable time in parametric, differential equations to demonstrate broader disease spread patterns . This method has been used to examine a wide range of phenomena—spanning all topics from the spread of measles and dengue fever to the efficacy of pulse vaccination policies [11, 12, 13]—but it cannot reveal the spread dynamics of any disease for which its key variables are unknown .
In contrast, agent-based models approach modeling from the level of individual organisms (agents) in heterogeneous environments, using the unique behaviors of individuals and the intra- and inter-species interactions which arise from these behaviors as modeling inputs [14, 15]. This input method allows for patterned behaviors of the whole complex system being modeled to emerge from the individual interactions over time. Using agent-based modeling, it is possible to capture real-world interactions between persons and other organisms more time-efficiently than is possible in the field and more immediately than is currently possible in the case of ZIKV with conventional differential equation models.
In this work, we developed an agent-based model to assess the impact of factors affecting the spread of ZIKV. The model observed ZIKV spread exclusively between humans and Ae. Aegypti mosquitoes, excluding the potential for person-to-person ZIKV spread since the risk chance of sexual transmission remains unknown . It provided evidence for a causal relationship between Zika outbreak severity and human population, leading us to hypothesize that Zika will, in the real world, most greatly afflict urban centers (as opposed to more rural areas) during locally transmitted outbreaks.
MATERIALS AND METHODS.
Our ZIKV model was implemented using Netlogo (version 5.3.1), the programming language and agent-based modeling environment created by U. Wilensky and colleagues [14, 17].
While creating the model studied in this work, we prototyped, incrementally developed, tested, and analyzed the model while also developing a parameter validation strategy unique to agent-based modeling. The model overall was validated by devising and valuing each parameter directly from the literature. Such literature focused mostly on the detailed behaviors of Aedes aegypti under variable environmental conditions and on ZIKV behavior in culture/laboratory settings.
Experiments were conducted approximating the environments of Rio de Janeiro State, Brazil and Wynwood, Florida in order to prove that the model could make predictions about ZIKV spread dynamics in locations of existing outbreaks [2, 3, 18, 19]. A set of simulations was also conducted and analyzed in two contexts to determine the factor which most greatly affected ZIKV outbreak severity in human agent populations. The following subsections outline experimental designs, parameter validations, and statistical analysis methods.
Modeled mosquitoes fed exclusively on human agents, excluding any possible divergence from a diet consisting solely of human blood [1, 20, 21, 22, 23]. Furthermore, only Ae. Aegypti were included while Ae. Albopictus were excluded since the former has demonstrated a higher capacity for ZIKV transmission and a tolerance for wider temperature ranges than has the latter [20, 21, 22, 23, 24]. Male Ae. Aegypti specimens were excluded and reproduction was abstracted in order to strictly restrict the diet of mosquito agents to human blood. Mosquito feeding was regulated by an energy variable which encouraged blood-feeding every 3-5 days . Initial energy values were stochastically determined from a range of 3-5 for each individual mosquito; this value was decremented daily starting 10 days after hatching (hatchlings remained fed until day 11) and decrementation continued until energy reached 0 (causing mosquito death) . Upon feeding, mosquitoes gained 2 units of energy, exchanged pathogens with its victim, and chanced random death (there was a 10% chance that the human agent would kill the mosquito when it landed; it was assumed that though the mosquito died in this case, pathogens were still exchanged). Maximum mosquito lifespans were stochastically, individually determined upon hatching from a range of 40-60 days; this range was based on Floridian average temperatures .
Mosquitoes could move a maximum of 20 meters per day starting on day 11 of their lives [26, 27]. Mosquitoes could sense and pursue human agents within this same 20 meter range; this chasing behavior was activated once the mosquitoes’ energy levels equaled 1-3. On days when the chase behavior was inactive or no human agents were in range, mosquitoes moved their maximum distance in a stochastically determined direction. Mosquitoes could not sense human agents after they were bitten greater than 20 times in a single day; this measure was implemented to control situations where humans were bitten thousands of times per day.
Humans could move, sense mosquito agents, and avoid them in a 4-meter radius. When avoiding mosquitoes, humans attempted to stochastically select a mosquito-free environment patch to move to. When no such a patch was available or when no mosquitoes were in range, humans moved their maximum range in a stochastically selected direction.
Mosquitoes or humans that acquired ZIKV did not become infectious until an extrinsic incubation period had passed. For humans, this period was 7 days; for mosquitoes, it was 10 days [1, 28, 29, 30; 31, 32, 33]. Since neither mosquitoes nor humans were classified as infectious until the concentration of ZIKV in their systems was high enough to have a very high chance of infecting the opposite species, there was a 100% chance of transmission when infected mosquitoes fed on uninfected humans and when uninfected mosquitoes fed on infected humans [1, 28, 29, 30, 31, 32, 33, 34].
Human agents lost their “ZIKV infectious” status 14 days after being initially bitten by a ZIKV infectious mosquito (7 days after humans gained the “infectious” status) . Their loss of infectious status was non-permanent; the possibility for acquired immunity is unsubstantiated in the literature and was therefore excluded. However, there is no evidence at the time of this writing to suggest that mosquitoes infected with ZIKV ever lose that infection; the “ZIKV infectious” designation was therefore irreversible in mosquito agents.
All simulated mosquito agents were female, but it was assumed that enough males were present to ensure successful egg fertilization, allowing simulated females to reproduce. Females had to be at least 11 days of age and be atop a standing water patch (the ideal egg-laying environment) [21, 26]. Since it was assumed that each female became engorged enough to reproduce after each feeding, they were allowed to reproduce (energy-wise) after each blood-meal . Females could reproduce 5 times total in their lives and laid 75 female (and therefore relevant) eggs per clutch [21, 26].
The starting spatial distribution of all human and mosquito agents was determined stochastically, placing all agents within the 1 km2 simulation environment. Standing water coverage was calculated as the ratio of standing water pixels to total pixels; these ratios were then converted to percentages. Collections of standing water were uniformly distributed throughout the simulation environment.
Experiment 1 studied the relationship between standing water coverage, the population density of humans (persons/km2), and ZIKV outbreak severity (as measured by the peak percentage of infected humans on any given day during the 100 day simulations). Standing water coverage incremented by 2% per simulation ranging from 2% to 20% while human population density incremented by 50 persons/km2 per simulation ranging from 179 to 529 persons/km2; each combination of standing water coverage and human population density within these intervals represented one simulation (or square in the figures) [Fig. 1, 2]. The initial mosquito population was arbitrarily set at 525 individuals/km2; the initial population of ZIKV infectious mosquitoes was set at 250 individuals/km2. All other variables remained constant at the values found in the literature. The dependent variable in experiment 1 was ZIKV outbreak severity.
Experiment 2 observed correlations between standing water coverage, human population density, and peak mosquito populations (as measured by the peak count of mosquito agents present on any given day during the 100 day simulations multiplied by 2 to account for male mosquitoes). Experiment 2 represents the same simulation set as experiment 1; the experiments merely draw from 2 different data outputs produced by the same simulation set. Therefore, the initial population of ZIKV infectious and non-ZIKV infectious mosquitoes in experiment 2 was the same as in experiment 1. The dependent variable in experiment 2 was peak mosquito population.
Experiment 3 sought to study ZIKV under approximations of conditions in Rio de Janeiro State, Brazil and the Wynwood neighborhood of Miami-Dade County, Florida. Two sets of simulations separate from the prior experiments were conducted: Florida’s set contained 16 trials while Brazil’s set contained 40 trials. The Florida simulations set standing water coverage equal to 22% and human population density equal to 2983 persons/km2 . The Brazil simulations set standing water coverage equal to 18% and human population density equal to 1174 persons/km2 [36, 37].
All experiments were conducted with simulations run on a 64-bit Windows 10 operating system with 16 GB of installed RAM. Experiments 1 and 2 are comprised of data from one set of 80 total simulations; each experiment explores the same data set in a different light. This 80-simulation set took approximately 72 hours to complete. Each simulation took on average 54 minutes to complete, but this average is computed cognizant that simulation lengths increased exponentially with human population density and standing water coverage. This is because both of these variables increased peak mosquito populations, which in turn directly increased processing demands. The data in Experiment 3 were obtained through 56 total simulations, amounting to roughly 50 hours of processing.
In experiment 1, a two-factor ANOVA test without replication was conducted to ascertain the relationship between human population density, standing water coverage, and outbreak severity [Table 1]. The null hypothesis for the rows, i.e. human population densities, was given by H0: there is no significant difference in Zika outbreak severity between the means of the human population. The null hypothesis for the columns, i.e. standing water coverages, was given by H0: there is no significant difference in Zika severity between the percentages of standing water coverage.
In experiment 2, a two-factor ANOVA test without replication was conducted to learn the relationship between human population density, standing water coverage, and the peak population of mosquitoes [Table 2]. The null hypothesis for the rows, i.e. human population densities, was given by H0: there is no significant difference in the peak mosquito population between the means of the human population. The null hypothesis for the columns, i.e. standing water coverages, was given by H0: there is no significant difference in the peak mosquito population between the percentages of standing water coverage.
Figure 1. Heat Map of ZIKV Outbreak Severities
Figure 1. Zika outbreak severity, given by the peak percentage of the human population infected with Zika in the first 100 days of simulation, is plotted in a heat map against standing water coverage percentages (x-axis) and human population densities (y-axis; humans/km2). Each colored square represents one simulation from the 80-simulation set whose human population density and standing water coverage percentage value settings are delineated by the x- and y-axes; the color gradient illustrates Zika outbreak severity.
Table 1. Results of ANOVA on Zika Outbreak Severity
We performed two-factor ANOVA test without replication on the simulated data to observe the relationship between human population density (rows), standing water coverage (columns), and Zika outbreak severity. The analysis results rejected the null hypothesis that there was no significant difference in Zika outbreak severity between the means of the human population. However, the analysis observed no significant difference in Zika outbreak severity between the percentages of standing water coverage.
Figure 1 plots Zika outbreak severity in each of the individual simulations within the experiment’s 80-simulation set on a heat map. This heat map declares human population density and standing water coverage percentage, both variables fixed at the start of each simulation, as the experiment’s independent variables. It asserts a causal relationship between these variables and Zika outbreak severity.
The ANOVA test in Table 1 further delineates this causal relationship. The test in Table 1 considered the relationships between simulated human population densities (rows), standing water coverage percentages (columns), and Zika outbreak severity [Table 1]. It revealed a significant difference in Zika outbreak severity between the means of the human population (p < 0.01) [Table 1]. It indicated no significant difference in Zika outbreak severity between the standing water coverages (p > 0.95) [Table 1].
Figure 2. Heat Map of Peak Mosquito Populations
Observing the relationship between the initial population density of humans (blood-meal availability; humans/km2), the amount of standing water present in the environment (breeding ground availability), and the peak population of simulated mosquitoes (given in units of 1,000 mosquitoes) via heat map. Each colored square represents one simulation from the same 80-simulation set as Experiment 1 whose human population density and standing water coverage percentage value settings are delineated by this figure’s x- (standing water coverage) and y- (human population density) axes; the color gradient illustrates peak mosquito populations.
Table 2. Results of ANOVA on Peak Mosquito Population
We also analyzed the relationship between human population density (rows), standing water coverage (columns), and the peak populations of mosquitoes with two-factor ANOVA. Both null hypotheses were rejected (see green highlighted p-values): there was a significant difference in peak mosquito populations between both the means of the human population and the percentages of standing water coverage.
Figure 2 visually illustrates the relationship between human population densities, standing water coverages, and peak mosquito populations in each of the simulations in the 80-simulation set [Fig. 2]. It asserts a causal relationship between the first two of these factors, the independent variables; and the third, the dependent variable. It claims that both standing water coverages (essentially the availability of mosquito breeding grounds) and human population densities (basically the availability of food for mosquitoes) affect the peak populations of mosquitoes.
The ANOVA test in Table 2 confirms the assertions put forth in Figure 2. The test in Table 2 looked at the relationships between human population densities (rows), standing water coverage percentages (columns), and peak mosquito populations [Table 2]. It confirmed a significant difference in the peak populations of mosquitoes between both the means of the human population and the percentages of standing water coverage (p < 0.01) [Table 2]. However, no significant differences between the peak populations of mosquitoes and outbreak severities were observed.
Figure 3. Outbreak Severity Comparison: Brazil vs. Florida
Simulated Zika outbreak severity under settings approximating the Wynwood neighborhood in Miami-Dade county, Florida (in orange) and Rio de Janeiro State, Brazil (in blue).
On average, Zika outbreak severity under settings approximating the Wynwood neighborhood in Miami-Dade county, Florida equaled 30.940% [Fig. 3]. Given that in Brazil, this average was 27.425%, Brazil conditions produced lower outbreak severities than Florida conditions. Since coefficient of variation corresponding with the average Florida outbreak severity was 6.471% (0.06471) while the that same coefficient corresponding with the average outbreak severity in Brazil was 10.826% (0.10826), more confidence can be placed in results from the Florida simulations than in results from the Brazil trials.
It is noteworthy that there was a significant difference (p < 0.01) between the simulated outbreak severities of Rio de Janeiro State and Wynwood. This difference indicates that in the absence of the disease control regulations in place in Florida, natural conditions in the Miami-Dade area may lend themselves towards worser outbreaks than Rio de Janeiro.
This study found a significant correlation (p < 0.01) between Zika outbreak severity and human population density in its simulations [Table 1]. We hypothesize that the significant correlation observed between Zika outbreak severity and human population density in this study translates into the real world of Zika spread. Therefore, these results are impactful on the broader public health field.
In future studies, we suggest further building on the model’s modular structure to make human agent behaviors more complex. Additionally, more minute deposits of standing water than were possible to account for on this model’s scale should be accounted for in future iterations. Seasonal and weather-related temperature and water deposit fluctuations should also be incorporated in future iterations. These additions and their added complexities should be accommodated with higher-powered hardware for running simulations.
When we began this study, we hypothesized that there would be a direct correlation between peak mosquito populations and outbreak severities; the tests in this study prove no such correlation. We now hypothesize that the observed lack of correlation is caused by the fact that while peak mosquito populations increased, the initial number of ZIKV infectious mosquitoes remained constant across simulations; as peak mosquito populations increased and the ratio of ZIKV infected mosquitoes to ZIKV uninfected mosquitoes decreased, mosquito-human infection likelihood decreased. We hypothesize further that this phenomenon was exacerbated by ZIKV uninfected mosquitoes beating out ZIKV infected mosquitoes for food and other resources. These revised hypotheses should be studied after future revisions to the model’s modular structure.
The model we produced was a novelty because of the early-designed behavioral stages that draw on both the limited available knowledge on the patterned behavior of ZIKV on a global scale and the more understood behavior of the disease’s vectors on a smaller scale. The results it produced contributed to our knowledge of the practical behavior of ZIKV. Most importantly, this model, in its success, served as a proof of concept for constructing similar agent-based disease models in future years.
This project was supported by the Department of Teaching and Learning at Vanderbilt University in addition to being supported by the Institute for Software Integrated Systems at Vanderbilt University and the School for Science and Math at Vanderbilt. Drs. Douglas Clark, Angela Eeds, Stephanie Weeden-Wright and innumerable others volunteered their time and mentorship in bringing this project to fruition.
- L.R. Petersen et al, Zika virus. The New England Journal of Medicine, 374, 1552-1563 (2016).
- CDC, “All countries & territories with active Zika virus transmission,” (2016; http://www.cdc.gov/zika/geo/active-countries.html).
- CDC, “Zika cases reported in the united states,” (2016; https://www.cdc.gov/zika/intheus/maps-zika-us.html#active-florida).
- J. Cerbino-Neto et al, Clinical manifestations of Zika virus infection Rio de Janeiro Brazil 2015, Emerging Infectious Diseases, 22, 2015–2017 (2016).
- U.S. National Institute of Neurological Disorders and Stroke, “Guillain-Barre Syndrome Fact Sheet,” (2014; http://www.ninds.nih.gov/disorders/gbs/detail_gbs.htm).
- D. Musso et al, Rapid spread of emerging Zika virus in the Pacific area, Clinical Microbiology and Infection, 20, O595-O596 (2014).
- J. Mlakar et al, Zika virus associated with Microcephaly, The New England Journal of Medicine, 374, 951-958 (2016).
- W. K. de Oliveira et al, Increase in reported prevalence of microcephaly in infants born to women living in areas with confirmed Zika virus transmission during the first trimester of pregnancy – Brazil 2015, Morbidity and Mortality Weekly Report, 65, 242-247 (2016).
- NSF, “National Science Foundation issues call for Zika virus proposals,” (2016; http://www.nsf.gov/news/news_summ.jsp?cntn_id=137621&org=NSF&from=news).
- D. Smith and L. Moore, “The SIR model for spread of disease – the differential equation model,” (2016; http://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-the-differential-equation-model).
- L. Stone, B. Shulgin, and Z. Agur, Theoretical examination of the pulse vaccination policy in the SIR epidemic model, Mathematical and Computer Modelling, 31, 207-215, (2000).
- Z. Feng and J.X. Velasco-Hernandez, Competitive exclusion in a vector-host model for the dengue fever, Journal of Mathematical Biology, 35, 523-544 (1997).
- O. N. Bjørnstad, B. F. Finkenstädt, and B. T. Grenfell, Dynamics of measles epidemics: estimating scaling of transmission rates using a time series SIR model, Ecological Monographs, 72, 169–184 (2002).
- U. Wilensky and W. Rand, An introduction to Agent-Based Modeling: Modeling Natural, Social, and Engineered Complex Systems with Netlogo (MIT Press, Cambridge, MA, 2015).
- C. M. Macal and M.J. North, Tutorial on agent-based modeling and simulation, Journal of Simulation, 4, 151-162 (2010).
- E. E. Petersen et al, Update: Interim Guidance for Preconception Counseling and Prevention of Sexual Transmission of Zika Virus for Persons with Possible Zika Virus Exposure – United States, September 2016, Morbitidy and Mortality Weekly Report, 65, 1077-1081 (2016).
- U. Wilensky, “Netlogo Homepage,” (2016; https://ccl.northwestern.edu/netlogo/).
- D. Chang, “Four more local Zika cases in Wynwood – the first one in Pinellas County,” Miami Herald (2016; http://www.miamiherald.com/news/health-care/article97345352.html).
- PAHO, “Zika epidemiological report: Brazil,” (2016; http://www.paho.org/hq/index.php?option=com_docman&task=doc_view&gid=35221<emid=270&lang=en).
- M. U. G. Kraemer et al, The global distribution of the arbovirus vectors Aedes aegypti and Aedes albopictus, Elife, 4, 1-18 (2015).
- C. Zettel and P. Kaufman, “Yellow fever mosquito – Aedes aegypti,” (Univ. of Florida, Gainesville, FL, 2016; http://entnemdept.ufl.edu/creatures/aquatic/aedes_aegypti.htm).
- CDC, “Zika Virus: Transmission & Risks,” (2016; www.cdc.gov/zika/transmission/).
- WHO, “Zika virus vectors and risk of spread in the WHO European Region,” (Geneva, Switzerland, 2016; http://www.euro.who.int/__data/assets/pdf_file/0007/304459/WEB-news_competence-of-Aedes-aegypti-and-albopictus-vector-species.pdf?ua=1).
- O. J. Brady et al, Modelling adult Aedes aegypti and Aedes albopictus survival at different temperatures in laboratory and field settings, Parasites & Vectors, 6, 351 (2013).
- G. A. H. McClelland and G.R. Conway, Frequency of blood feeding in the mosquito Aedes aegypti, Nature, 232, 485-486 (1971).
- CDC, “Dengue: mosquito life cycle,” (2016; http://www.cdc.gov/dengue/resources/factSheets/MosquitoLifecycleFINAL.pdf).
- P. T. McDonald, Population characteristics of domestic Aedes aegypti (diptera: culicidae) in villages on the Kenya coast ii. dispersal within and between villages, Journal of Medical Entomology, 14, 49–53 (1977).
- European Centre for Disease Prevention and Control, “Factsheet for health professionals: Zika Virus Disease,” (Health Topics, 2016; http://ecdc.europa.eu/en/healthtopics/zika_virus_infection/factsheet-health-professionals/Pages/factsheet_health_professionals.aspx).
- NZ Ministry of Health, About Zika Virus Disease, Diseases and Conditions (2016; http://www.health.govt.nz/our-work/diseases-and-conditions/zika-virus).
- WHO, “Zika Virus Fact Sheet,” WPRO Fact Sheets (2016; http://www.wpro.who.int/mediacentre/factsheets/fs_05182015_zika/en/).
- H.L.C. Dutra et al, Wolbachia blocks currently circulating Zika virus isolates in Brazilian Aedes aegypti mosquitoes, Cell Host & Microbe, 19, 1-4 (2016).
- M. Dubrulle et al, Chikungunya virus and Aedes Mosquitoes: saliva is infectious as soon as two days after oral infection, PLoS ONE, 4, e5895 (2009).
- M. I. Li et al, Oral susceptibility of singapore Aedes (Stegomyia) aegypti (Linnaeus) to Zika virus, PLoS Neglected Tropical Diseases, 6, (2012).
- A.M. Oster et al, Update: interim guidance for prevention of sexual transmission of Zika virus – United States 2016, Morbidity and Mortality Weekly Report, 65, 323-325 (2016).
- U.S. Census Bureau, “2015 U.S. Gazetteer Files,” (2015; http://www.census.gov/geo/maps-data/data/gazetteer2015.html).
- R. M. Schneider et al, “Rio de Janeiro Brazil,” (Encyclopedia Britannica, 2016; https://www.britannica.com/place/Rio-de-Janeiro-Brazil).
- U.N. Data, “City population by sex, city and city type,” (2016; http://data.un.org/Data.aspx?d=POP&f=tableCode%3A240).