# Epidemiological significance of dengue virus genetic variation in mosquito infection dynamics

Albin Fontaine, Sebastian Lequime, Isabelle Moltini-Conclois, Davy Jiolle, Isabelle Leparc-Goffart, Robert Charles Reiner, Louis Lambrechts

► To cite this version:

Albin Fontaine, Sebastian Lequime, Isabelle Moltini-Conclois, Davy Jiolle, Isabelle Leparc-Goffart, et al.. Epidemiological significance of dengue virus genetic variation in mosquito infection dynamics. PLoS Pathogens, 2018, 14 (7), pp.e1007187. 10.1371/journal.ppat.1007187. hal-01858750

HAL Id: hal-01858750

https://hal.science/hal-01858750

Submitted on 22 Aug 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Epidemiological significance of dengue virus genetic variation in mosquito infection dynamics

Albin Fontaine¹,²,³,⁴*, Sebastian Lequime¹,³,⁵, Isabelle Moltini-Conclois¹,³, Davy Jiolle¹,³, Isabelle Leparc-Goffart⁶, Robert Charles Reiner, Jr⁷, Louis Lambrechts¹,³*

1 Insect-Virus Interactions Group, Department of Genomes and Genetics, Institut Pasteur, Paris, France,

2 Unité de Parasitologie et Entomologie, Département des Maladies Infectieuses, Institut de Recherche Biomédicale des Armées, Marseille, France, 3 Génomique Evolutive, Modélisation et Santé, Unité Mixte de Recherche 2000, Centre National de la Recherche Scientifique, Paris, France, 4 Aix Marseille Université, IRD, AP-HM, SSA, UMR Vecteurs-Infections Tropicales et Méditerranéennes (VITROME), IHU-Méditerranée Infection, Marseille, France, 5 KU Leuven Department of Microbiology and Immunology, Rega Institute, Laboratory of Clinical and Epidemiological Virology, Leuven, Belgium, 6 Centre National de Référence des Arbovirus, Institut de Recherche Biomédicale des Armées, Marseille, France, 7 Institute for Health Metrics and Evaluation, University of Washington, Seattle, Washington, United States of America

* albinfont@gmail.com (AF); louis.lambrechts@pasteur.fr (LL)

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

## OPEN ACCESS

**Citation:** Fontaine A, Lequime S, Moltini-Conclois I, Jiolle D, Leparc-Goffart I, Reiner RC, Jr, et al. (2018) Epidemiological significance of dengue virus genetic variation in mosquito infection dynamics. PLoS Pathog 14(7): e1007187. https://doi.org/10.1371/journal.ppat.1007187

**Editor:** Neil M. Ferguson, Imperial College London, UNITED KINGDOM

**Received:** February 14, 2018

**Accepted:** June 28, 2018

**Published:** July 13, 2018

**Copyright:** © 2018 Fontaine et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** This work was supported by Agence Nationale de la Recherche (grants ANR-16-CE35-0004-01 and ANR-17-ERC2-0016-01), the French Government's Investissement d'Avenir program Laboratoire d'Excellence Integrative Biology of Emerging Infectious Diseases (grant ANR-10-LABX-62-IBEID), the City of Paris Emergence(s) program in Biomedical Research and U.S. National

## Abstract

The kinetics of arthropod-borne virus (arbovirus) transmission by their vectors have long been recognized as a powerful determinant of arbovirus epidemiology. The time interval between virus acquisition and transmission by the vector, termed extrinsic incubation period (EIP), combines with vector mortality rate and vector competence to determine the proportion of infected vectors that eventually become infectious. However, the dynamic nature of this process, and the amount of natural variation in transmission kinetics among arbovirus strains, are poorly documented empirically and are rarely considered in epidemiological models. Here, we combine newly generated empirical measurements *in vivo* and outbreak simulations *in silico* to assess the epidemiological significance of genetic variation in dengue virus (DENV) transmission kinetics by *Aedes aegypti* mosquitoes. We found significant variation in the dynamics of systemic mosquito infection, a proxy for EIP, among eight field-derived DENV isolates representing the worldwide diversity of recently circulating type 1 strains. Using a stochastic agent-based model to compute time-dependent individual transmission probabilities, we predict that the observed variation in systemic mosquito infection kinetics may drive significant differences in the probability of dengue outbreak and the number of human infections. Our results demonstrate that infection dynamics in mosquitoes vary among wild-type DENV isolates and that this variation potentially affects the risk and magnitude of dengue outbreaks. Our quantitative assessment of DENV genetic variation in transmission kinetics contributes to improve our understanding of heterogeneities in arbovirus epidemiological dynamics.
Institutes of Health grant 1P01AI098670-01A1. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

## Author summary

Transmission of arboviruses by mosquitoes is a time-dependent process because dissemi-nation of the virus from the insect's digestive tract, where the infection is initially estab-lished, to the salivary glands, where it can be transmitted to the next host, typically takes several days. The duration of this time interval is an important determinant of the propor-tion of infected mosquitoes that live long enough to become infectious. Due to the lack of available data, epidemiological models rarely account for the dynamic nature of this pro-cess and its natural variation among arbovirus strains. In this study, we experimentally monitored the dynamics of mosquito infection by eight strains of dengue virus and found significant differences between them. In addition, we used the empirical data obtained in a simulation model to show that the differences observed among the strains are expected to result in human dengue outbreaks of different probability and size. Together, our results indicate that natural variation in mosquito infection dynamics between arbovirus strains may contribute to the unexplained heterogeneity of dengue virus transmission patterns.

## Introduction

In the last few decades, arthropod-borne viruses (arboviruses) have become major contribu-tors of global mortality and disability, calling for an improved understanding of their transmis-sion dynamics [1]. Transmission of arboviruses is governed by a complex set of biotic and abiotic factors [2], the majority of which relate to the vector biology [3]. Individual, temporal and spatial variability in these factors results in considerable heterogeneities in arbovirus trans-mission dynamics, but the relative contribution of each factor is still poorly understood [4-7].

In classic models of mosquito-borne pathogen transmission, one of the most influential parameters is the extrinsic incubation period (EIP), defined as the time interval between path-ogen acquisition and pathogen transmission by the vector [8]. In these models, EIP interacts with the probability of mosquito survival as an exponential term, so that even a small deviation in EIP will substantially change the number of mosquitoes surviving long enough to become infectious. Despite the recognized importance of EIP as a powerful determinant of mosquito-borne pathogen transmission, variation in the duration of EIP has rarely been accounted for in models of mosquito-borne pathogen epidemiology [9]. Current approaches to assess the epi-demiological significance of variation in arbovirus transmission kinetics rely on temperature-dependent models developed from limited empirical data [10-12]. These models generally assume that the only factor affecting transmission kinetics is ambient temperature.

Here, we investigated whether arbovirus genetic variation may alter transmission kinetics in epidemiologically significant ways. Specifically, we experimentally measured variation in transmission kinetics among several wild-type isolates of dengue virus (DENV) type 1 repre-senting the worldwide genetic diversity of recently circulating strains. DENV are positive-sense, single-stranded RNA viruses in the *Flavivirus* genus that cause more human infections than any other arbovirus [13]. They exist as four genetically divergent types that loosely cluster antigenically [14]. Each DENV type is divided into major phylogenetic clades, often termed genotypes, which in turn contain a wide diversity of strains [15]. All four serotypes of DENV are mainly transmitted by the mosquito *Aedes aegypti* worldwide [16].

A mosquito's intrinsic ability to transmit a virus, referred to as vector competence, is a step-wise procedure that begins when the mosquito takes an infectious blood meal on a viremic host. Following initial infection and replication in the midgut, the virus disseminates systemically from the midgut to secondary tissues including salivary glands, where it can be released
into saliva during a subsequent blood meal. The probability of detecting DENV in *Ae. aegypti* saliva is correlated with the amount of virus that has spread out of the midgut [17], which is often used as a proxy for transmission potential. Although it is generally measured by end-point laboratory assays, vector competence is a dynamic process that combines with vector survival and EIP to determine the proportion of infectious vectors over time [18]. The duration of DENV EIP is in the range of 8–12 days at a temperature of 25–28°C [12, 19, 20], which is close to adult mosquito life expectancy under field conditions.

The influence of EIP on DENV epidemiology was recently examined in the context of *Wolbachia*-mediated blocking of virus transmission [21, 22]. Earlier studies suggested that viral genetic variation in EIP could also influence DENV epidemiology in nature. For example, differences of EIP among DENV type 2 strains have been proposed as a mechanism of genotype replacement and dengue emergence [23]. By comparing three DENV type 2 strains of the American genotype with three DENV type 2 strain of the Southeast Asian genotype, Anderson and Rico-Hesse estimated a 2- to 65-fold increase in vectorial capacity for the invasive Southeast Asian genotype relative to the resident American genotype [23]. This increase resulted from a 7-day difference in EIP between genotypes (7 days for the Southeast Asian vs. 14 days for the American genotype) across a range of daily survival rates. However, EIP was estimated as a single ‘threshold’ value defined as the time at which 50% of salivary glands became infected, which likely misrepresents EIP for an entire population [24].

To more accurately assess the epidemiological significance of DENV genetic variation in transmission kinetics, we used two complementary approaches. First, we experimentally monitored the cumulative proportion of mosquitoes with a systemic infection over time for eight genetically diverse, low-passage DENV type 1 isolates. We then used these *in vivo* empirical data to derive a 3-parameter model of the kinetics of systemic infection, a proxy for transmission potential, and compare transmission dynamics among DENV isolates. Second, we implemented a stochastic agent-based model accounting for intra-mosquito and intra-human viral dynamics to compute time-dependent individual transmission probabilities and simulate arbovirus outbreaks *in silico*. Using this model, we showed that the empirically observed differences in DENV transmission kinetics have the potential to significantly influence the probability and size of human dengue outbreaks.

## Results

### Diverse field-derived DENV isolates display variation in systemic infection kinetics

We orally exposed wild-type *Ae. aegypti* mosquitoes to a standardized infectious dose of eight DENV type 1 isolates derived from human sera sampled in various endemic areas worldwide during 2010–2013 (Table 1). The eight DENV type 1 isolates were genetically diverse and included six isolates belonging to genotype I, and two isolates belonging to genotype V (Fig 1). Mosquitoes were collected 4, 6, 8, 12 and 18 days after virus oral exposure and assayed for DENV midgut and/or systemic infection. For each isolate, the mean number of mosquitoes assayed per time point was 25 and ranged from 10 to 32.

We first tested infection status (S1 Fig) and subsequently focused on systemic infection dynamics of midgut-infected mosquitoes. Systemic DENV infection was determined at each time point based on virus presence in head tissues. Systemic infection prevalence among midgut-infected mosquitoes ranged from 4.5% to 100% across isolates and time points and was significantly influenced by the time point (p = 2.0×10⁻¹⁶) but not by uncontrolled variation in blood meal titers (p = 0.08) or the interaction between time point and blood meal titer (p = 0.61).
<table><caption>**Table 1.** Description of DENV type 1 isolates used in this study. The blood meal titer refers to the concentration of infectious viral particles (expressed in log₁₀-transformed focus-forming units per mL) measured in the artificial blood meal offered to mosquitoes. The passage history refers to the number of prior amplifications in C6/36 cells.</caption><thead><tr><td>DENV-1 isolate</td><td>Genotype</td><td>Geographical origin</td><td>Year of isolation</td><td>Blood meal titer (log₁₀ FFU/mL)</td><td>Passage history</td><td>GenBank accession number</td></tr></thead><tbody><tr><td>Thai2010a</td><td>I</td><td>Thailand</td><td>2010</td><td>5.74</td><td>5</td><td>HG316481</td></tr><tr><td>Thai2010b</td><td>I</td><td>Thailand</td><td>2010</td><td>5.70</td><td>6</td><td>HG316482</td></tr><tr><td>Thai2012</td><td>I</td><td>Thailand</td><td>2012</td><td>5.80</td><td>3</td><td>MG877554</td></tr><tr><td>Thai2013</td><td>I</td><td>Thailand</td><td>2013</td><td>5.80</td><td>2</td><td>MG877556</td></tr><tr><td>Laos2012</td><td>I</td><td>Laos</td><td>2012</td><td>5.85</td><td>3</td><td>MG877552</td></tr><tr><td>NCal2013</td><td>I</td><td>New Caledonia</td><td>2013</td><td>5.77</td><td>2</td><td>MG877555</td></tr><tr><td>Gabon2012</td><td>V</td><td>Gabon</td><td>2012</td><td>5.82</td><td>5</td><td>MG877557</td></tr><tr><td>Haiti2012</td><td>V</td><td>Haiti</td><td>2012</td><td>5.81</td><td>3</td><td>MG877553</td></tr></tbody></table>

https://doi.org/10.1371/journal.ppat.1007187.t001

To quantify the kinetics of systemic DENV infection, we assumed that the cumulative change in the proportion of mosquitoes with a systemic DENV infection over time post exposure followed a sigmoid function. A logistic model (see *Materials and Methods*) was fitted to the empirical data to estimate three parameters (*K*, *B*, *M*) describing the kinetics of each isolate (Table 2; Fig 2A). *K* is the saturation level and represents the maximum proportion of infected mosquitoes with a systemic infection. *B* is the slope factor and represents the maximum slope of the cumulative function scaled by *K*. *M* is the lag time and represents the time at which the absolute increase in cumulative proportion of systemic infection is maximal, which is also when the systemic infection prevalence equals *K*/2.

Among isolates, the saturation level *K* ranged from 79% to 100% (mean ± SD: 93 ± 7%). The slope factor *B* was transformed into Δ*t*, the time required to rise from 10% to 90% of the saturation level, which ranged from 1.85 to 7.73 days (mean ± SD: 4.37 ± 2.03 days). The lag time *M* ranged from 4.84 to 7.36 days (mean ± SD: 5.68 ± 0.82 days). There was no statistically significant correlation between the three parameters (S2 Fig). However, isolates with a high saturation level tended to display a large slope factor and a short lag time. To graphically represent EIP variation among DENV isolates, we used the distribution of lag time (*M*) values as a proxy. We obtained the daily rate of new systemic infections by calculating the first derivative of the cumulative systemic infection function (Fig 2B). The curve of this rate over time is equivalent to the distribution of lag time values among individual mosquitoes.

To quantitatively compare the kinetics of systemic infection among DENV isolates, we fitted a global likelihood model based on the 3-parameter logistic function to the empirical data for all combinations of permutations of isolates (see *Materials and Methods*). Based on Akaike Information Criterion (AIC) values, we identified four groups of isolates (Fig 3B and 3C). The first group had slow systemic infection dynamics and included three DENV isolates (Gabon2012, NCal2013, Thai2012) with a relatively low saturation level, small slope factor and long lag time. The second group had fast systemic infection dynamics and included three DENV isolates (Thai2010b, Thai2013, Thai2010a) with a higher saturation level, larger slope factor and shorter lag time. The last two groups consisted of one DENV isolate each (Laos2012 and Haiti2012) that had unique combinations of parameters (Table 2). A hierarchical clustering analysis performed on the affinity matrix derived from AIC scores of all combinations of isolates revealed that isolate Laos2012 branched with the fast-disseminating infection group whereas isolate Haiti2012 branched with the slow-disseminating infection group (Fig 3C). The groups did not relate to DENV type 1 genotypes, as for example Haiti2012 and Gabon2012 (both from genotype V) did not group together.

The hierarchical clustering results for midgut infection prevalence did not match the results obtained for systemic infection prevalence, indicating that midgut infection patterns are poor
<figure><img src="image_4.png" /><figcaption>**Fig 1. Phylogenetic relationships between DENV isolates of the study.** Full-length genome sequences were used to infer the phylogenetic relationships among the isolates of the study in a background of representative sequences of DENV type 1. Genotypes are represented with distinct colors and the eight DENV isolates of the study are shown in black. Inferences were calculated with a maximum-likelihood method implemented in RAxML v.8.2.10[50]. The scale bar represents 0.007 substitutions per site.</figcaption></figure>

https://doi.org/10.1371/journal.ppat.1007187.g001

predictors of differences in the kinetics of systemic infection. Likewise, there was no association between the viral phylogeny and the kinetics of systemic infection (S3 Fig). As expected with such a small number of isolates, there was no statistically significant association between individual DNA polymorphisms or amino-acid differences with any of the three parameters describing the kinetics of systemic infection.
<table><caption>**Table 2.** Parameter estimates of systemic infection kinetics for each DENV isolate. Parameter values for each isolate were inferred from the cumulative change in the proportion of mosquitoes with a systemic DENV infection over time post exposure using model optimization. *K* is the saturation level and represents the maximum proportion of mosquitoes with a systemic infection. *B* is the slope factor and represents the maximum slope of the cumulative function scaled by *K*. Δ*t* is derived from *B* and represents the time required to rise from 10% to 90% of the saturation level. *M* is the lag time and represents the time at which the absolute increase in cumulative proportion is maximal, which is also when the systemic infection prevalence equals *K*/2. Background shading (grey scale) indicates the best isolate groupings as determined by the AIC method.</caption><thead><tr><th>Isolate</th><th><i>K</i> (%)</th><th><i>B</i> (% day<sup>-1</sup>)</th><th>Δ<i>t</i> (days)</th><th><i>M</i> (days)</th></tr></thead><tbody><tr><td>Gabon2012</td><td>93</td><td>68</td><td>6.45</td><td>5.92</td></tr><tr><td>NCal2013</td><td>87</td><td>107</td><td>4.11</td><td>5.96</td></tr><tr><td>Thai2012</td><td>79</td><td>57</td><td>7.73</td><td>5.37</td></tr><tr><td>Haiti2012</td><td>94</td><td>92</td><td>4.78</td><td>7.36</td></tr><tr><td>Thai2010b</td><td>98</td><td>124</td><td>3.55</td><td>5.15</td></tr><tr><td>Thai2013</td><td>100</td><td>96</td><td>4.57</td><td>5.92</td></tr><tr><td>Thai2010a</td><td>96</td><td>227</td><td>1.94</td><td>4.84</td></tr><tr><td>Laos2012</td><td>100</td><td>237</td><td>1.85</td><td>4.90</td></tr></tbody></table>

https://doi.org/10.1371/journal.ppat.1007187.t002

## Empirically observed variation in DENV transmission kinetics has the potential to dramatically impact the risk and magnitude of dengue outbreaks

To assess the epidemiological significance of the differences observed in the kinetics of DENV systemic infection in mosquitoes, we implemented a stochastic agent-based model that accounted for variation in the three parameters (*K*, *B*, *M*) of the logistic model described above, and their effect on the probability and size of dengue outbreaks. This method allows individual humans and mosquitoes, each one being characterized by a unique set of infection parameters, to interact within a simulated environment, and was designed to model the unfolding of real-world arboviral epidemics (see *Materials and Methods*). The 3-parameter function that was used to model the kinetics of systemic mosquito infection was incorporated as a time-dependent probability function to determine transmission status when a contact occurred between a human and a mosquito. Starting with one infected human in a population of 10,000 individuals, the model was run 100 independent times over 400 days using different combinations of the three parameters in the same range as the empirical values. The total number of infected humans was recorded for each combination of realistic parameter values. All combinations of parameter values resulted in a certain proportion of outbreaks without any secondary human infection (S4 File). Outbreaks with ≥1 secondary human infections were mainly partitioned into small-scale outbreaks (<100 secondary infections) and large-scale outbreaks (≥100 secondary infections) leading to a strongly bimodal distribution of outbreak size (S4A Fig). Lag time (*M*) had the most significant effect on the initiation of outbreaks (S4B and S4C Fig).

Based on the empirical values of their parameter combinations, the risk and magnitude of dengue epidemics varied among the eight DENV isolates (Fig 4A and 4B). The probability that at least 1 secondary human infection occurred ranged from 49% to 68%, and the probability of a large-scale dengue outbreak (≥100 secondary human infections) ranged from 38% to 50% among isolates (Fig 4A). During large-scale outbreaks, the mean proportion of the human population that became infected ranged from 69% to 90% among isolates (Fig 4B). In particular, the Haiti2012 isolate, which had the longest lag time (*M*) value, consistently led to a smaller outbreak size than all other isolates (Fig 4B). Overall, the probability and size of dengue outbreaks tended to be more similar among isolates from the same AIC-based group of infection kinetics described above. Differences between DENV isolates did not change significantly
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig 2.** Variation in the kinetics of systemic mosquito infection between DENV isolates. (A) The cumulative prevalence of systemic infection over time post virus exposure is shown for each DENV isolate. Data points represent the observed prevalence at each time point with their size being proportional to the sample size (number of mosquitoes tested). Dashes represent the 95% confidence interval of the prevalence. Lines correspond to the fitted values obtained with a 3-parameter logistic model. (B) The daily rate of new systemic infections (instantaneous velocity) over time post virus exposure is shown for each DENV isolate. It was calculated as the first derivative of the cumulative systemic infection function and is equivalent to the frequency distribution of lag time values among individual mosquitoes. In both graphs, the lower right panel shows the merged fitted values for all isolates.</figcaption></figure>

https://doi.org/10.1371/journal.ppat.1007187.g002
<figure><img src="image_7.png" /><figcaption>**Fig 3.** Comparison of systemic mosquito infection kinetics among DENV isolates. A logistic 3-parameter model was fitted to all permutations of isolates, with isolates from the same group being forced to share the same parameters. Panel (A) represents the isolate grouping that maximizes the probability function for systemic infection dynamics. This isolate grouping had the lowest AIC value among all tested permutations (S5 Fig). For each isolate group, data points represent the empirical measurements color-coded by isolate and the black line is the logistic fit of the group. An affinity matrix among isolates was derived from AIC values obtained from all model comparisons (S5 Fig). Panel (B) shows the dendrogram obtained from a hierarchical clustering analysis performed on the affinity scores. The distance between isolates corresponds to the number of times two isolates shared the same groupings with respect to their AIC. The red dotted line shows the cut-off value leading to the best AIC-based isolate grouping.</figcaption></figure>

https://doi.org/10.1371/journal.ppat.1007187.g003

when EIP was modeled by a continuous 3-parameter function in individual mosquitoes or by a threshold function with variable EIP among mosquitoes (S1 File). However, when a single EIP value was assigned to the entire mosquito population it resulted in significant discrepancies between the simulations outputs (S1 File).

## Discussion

Building on earlier studies that showed EIP variation between DENV genotypes [23] and possibly serotypes [19, 20, 25], we have demonstrated that viral genetic variation in the kinetics of mosquito infection may contribute to heterogeneity in DENV transmission dynamics. Such temporal variation in the transmission process has generally been ignored in models of mosquito-borne pathogen epidemiology [9], with the exception of temperature-dependent models [10–12]. We experimentally compared field-derived DENV isolates representing the worldwide genetic diversity of recently circulating strains. We found significant variation in the dynamics of systemic infection of *Ae. aegypti* among the DENV isolates, which was poorly predicted by midgut infection patterns. Using large-scale simulations with a stochastic agent-based model, we illustrated how the empirical variation in systemic infection dynamics may translate into sizeable differences in epidemiological outcome (risk and magnitude of dengue outbreak in the human population).

While previous studies established a link between empirical differences in EIP duration and epidemiological fitness of arbovirus strains [23, 26], the effect of EIP on transmission has
<figure><img src="image_8.png" /><figcaption>**Fig 4.** Simulated effect of variation in mosquito infection dynamics on the risk and magnitude of dengue outbreaks. A stochastic agent-based model was run 100 times for each combination of mean values of the three parameters (K, B, M) describing the kinetics of systemic mosquito infection. Other parameters, such as relative mosquito density, biting rate and human viral dynamics were held constant. Panel (A) shows the proportion of simulations that resulted in ≥100, <100 and no secondary human infections using the K, B and M values empirically measured for the eight DENV isolates. Panel (B) shows the total number of humans that became infected during large-scale dengue outbreaks (≥100 secondary human infections) using the empirical K, B and M values of the eight DENV isolates. In both panels, isolates are color-coded according to their AIC-based group of systemic mosquito infection kinetics.</figcaption></figure>

https://doi.org/10.1371/journal.ppat.1007187.g004

typically been extrapolated using a single value of EIP for the entire mosquito population. Our study went one step further by considering EIP as a quantitative trait following a population distribution [27]. Our results support the conclusion that the classic approach of measuring EIP as a 'threshold' time point at which the first mosquitoes become infectious is likely misleading for a mosquito population (S1 File). We used the cumulative proportion of infectious mosquitoes over time to infer the population distribution of lag time values (our proxy for EIP). Our results confirm that the distribution of lag time values cannot be described with a fixed threshold model because it derives from the sigmoid shape of a logistic function [28]. Furthermore, the inferred distribution of individual EIP values had a different shape between DENV isolates, adding another layer of complexity for comparing the epidemiological and evolutionary potential of different strains. Thus, our results lead to the conclusion that EIP is a quantitative trait that cannot be encapsulated in a single value for a given set of conditions. Although our general recommendation for future modeling studies is to use a full 3-parameter model whenever possible, EIP₅₀ can be considered as the most relevant single parameter because M (half of the maximum proportion of infected mosquitoes with a systemic infection) was the most influential in our model.

It is well established that the duration of DENV EIP can be influenced by environmental factors such as temperature [12, 20, 29]. Recently, Ye and colleagues used a breeding design to demonstrate that the duration of DENV EIP is also influenced by mosquito genes [30]. Genetic variation is a prerequisite for the evolution of a trait, but it is unlikely that natural selection would act on mosquito populations to drive evolution of DENV EIP. This is because DENV infection in natural *Ae. aegypti* populations is rare [31] and does not significantly affect mosquito lifespan [32]. Conversely, variation in EIP duration among DENV strains could fuel the evolution of a shorter EIP. Under field conditions, a significant proportion of adult mosquitoes are expected to die before they are capable of transmitting virus. A shorter EIP would confer virus variants a significant evolutionary advantage by increasing their probability of transmission. However, DENV transmission rarely occurs before at least 7 days of EIP in *Ae. aegypti*.
[12, 19, 20, 29]. Unlike other arboviruses with an EIP that can be as short as 2–3 days [33, 34], the longer duration of DENV EIP is presumably evolutionarily constrained. Results from the present study indicate that evolution of a shorter DENV EIP is not restricted by a lack of standing genetic variation for this trait among virus strains.

Despite the evidence for genetic variation among virus variants, we did not detect any association between viral nucleotide polymorphisms and the parameters describing the kinetics of systemic mosquito infection. There was no consistent association between kinetics and DENV genotypes. For example, a 2012 isolate from Thailand belonging to genotype I displayed the longest lag time and a slow-disseminating profile. This isolate was phylogenetically closest to a 2012 isolate from Laos within genotype I, which displayed a fast-disseminating profile. Small nucleotide differences in the genome sequence can disproportionately affect DENV replication by modifying the secondary structure of the viral RNA molecule [35], which may contribute to large phenotypic differences between closely related viruses. However, the most likely explanation for the lack of genotype-phenotype association is that our study was underpowered for a genetic association analysis due to the small number of isolates examined.

We used the kinetics of systemic mosquito infection as a proxy for EIP dynamics due to the logistical and ethical constraints of experimentally measuring mosquito-to-human DENV transmission. We assumed that the virus presence in mosquito head tissues indicated the potential for onward transmission based on evidence that DENV is detected a few days earlier in head tissues than in salivary glands [36], and lack of evidence for a salivary gland escape barrier in *Ae. aegypti* [37]. However, it is possible that our measurements of EIP duration based on the lag time of systemic infection are underestimates. The probability of detecting DENV in *Ae. aegypti* saliva is positively correlated with the infectious titer of disseminated virus [17]. Therefore, we cannot exclude that some mosquitoes had virus titers in their head tissues that were high enough to be detected, but too low to result in virus release in saliva. The extent of underestimation is unknown but should not change the interpretation of differences between strains under the assumption that this underestimation is proportionate. Another limitation of our study is that we used a very narrow range of infectious doses to orally challenge the mosquitoes. Infectious dose is one of the most powerful determinants of DENV infection probability and EIP duration in *Ae. aegypti* [22, 38, 39]. It is therefore possible that the parameter estimates of mosquito infection kinetics would have been different if we had used a larger range of infectious doses. Finally, we used a minimal number of cell culture passages to produce the viral stocks, which may have promoted adaptation to mosquito cells. However, we know of no evidence that adaptation to mosquito cell culture would translate into differences in viral dissemination *in vivo*.

One important aspect of our results is that we evaluated the epidemiological significance of the observed variation in transmission kinetics among DENV isolates. We used mathematical modeling to incorporate variation in the dynamics of DENV transmission by mosquitoes into a relevant epidemiological framework. We implemented an agent-based model that included stochasticity for different parameters, such as the daily survival rate of mosquitoes and the probability of human-mosquito contact. Importantly, our model considered intra-host infection dynamics, both in humans and mosquitoes. Consistent with results from deterministic Ross-MacDonald models of vector-borne pathogen transmission [8, 40], our simulations indicated that EIP duration was the most influential factor underlying DENV transmission dynamics. In particular, EIP duration had a stronger influence than the maximum proportion of infectious mosquitoes on both the probability and the size of dengue outbreaks. We emphasize that these differences in transmission dynamics depend on the validity of our modeling assumptions. For instance, there could be other differences between DENV strains, such as human viral dynamics, that could counteract or modulate the effects observed in our model.
In conclusion, this study illustrated the power of combining empirical data obtained from experiments *in vivo* with epidemiological modeling to unravel the determinants of vector-borne pathogen transmission dynamics. DENV transmission dynamics are characterized by heterogeneities resulting from individual, temporal and spatial variabilities that are often difficult to identify [4–7]. We have shown that variability in mosquito infection dynamics exists among arbovirus strains and is epidemiologically meaningful. Such quantitative dissection of transmission dynamics contributes to improve our understanding of heterogeneities in arbovirus epidemiology.

# Materials and methods

## Experimental approach

**Ethics statement.** The Institut Pasteur animal facility has received accreditation from the French Ministry of Agriculture to perform experiments on live animals in compliance with the French and European regulations on care and protection of laboratory animals. This study was approved by the Institutional Animal Care and Use Committee at Institut Pasteur under protocol number 2015–0032.

**Dengue virus isolates.** Eight DENV type 1 isolates derived from human sera were selected from the collection of the French National Reference Center for Arboviruses based on their geographical origin and year of isolation in order to represent the worldwide viral genetic diversity of recently circulating strains (Table 1). Informed consent of the patients was not necessary because virus isolated in laboratory cell culture is no longer considered a human sample. Prior to their use in mosquito experimental infections, DENV isolates were amplified 2–6 times in C6/36 (*Aedes albopictus*) cells. Virus stock was prepared as previously described [27]. Virus titration was performed by focus-forming assay (FFA) in C6/36 cells after 5 days of incubation at 28°C, as previously described [27].

**Experimental mosquito infections.** Experiments were carried out with wild-type *Ae. aegypti* mosquitoes derived from a field population originally sampled in the Muang District of Kamphaeng Phet Province, Thailand and took place within 10 generations of colonization. Mosquitoes were maintained and orally exposed to DENV as previously described [27]. Briefly, five- to seven-day-old females deprived of sucrose solution for 24 hours were allowed to feed for 15 min from a membrane feeding system (Hemotek) that contained an artificial infectious blood meal maintained at 37°C. The infectious blood meal consisted of two volumes of washed rabbit erythrocytes and one volume of pre-diluted viral stock, supplemented with 10 mM adenosine triphosphate as a phagostimulant. The viral stock was pre-diluted in order to standardize the infectious dose among the DENV isolates. The infectious dose was chosen to result in an expected infection rate of 50–100%. The same infectious blood meal was used to expose three sequential batches of mosquitoes during a total of 45 min. An aliquot of the artificial blood meal was collected immediately prior and after blood feeding and stored at -80°C for later titration by FFA. The mean (± SD) log₁₀-transformed viral titer dropped from 5.73 ± 0.13 to 5.46 ± 0.26 FFU/ml during the 45-min feeding period. This difference was statistically significant (t test: *p* = 0.0093) so the average of the two measurements for each isolate was used in all subsequent analyses. After virus exposure, fully engorged females were transferred to 1-pint cardboard containers and maintained under controlled conditions as described above with a maximum of 30 females per container.

**DENV detection.** The presence of a midgut and/or systemic DENV infection was determined on day 4, 6, 8, 12 and 18 post virus exposure. An infected body indicates a midgut infection, and an infected head indicates a systemic (disseminated) infection. Batches of mosquitoes were freeze-killed at each time point and their heads and bodies were separated and
stored at -80°C until processing. Virus detection was performed for all mosquito bodies and for the heads of all mosquitoes with a DENV-positive body. The majority of samples were tested using a qualitative version of the FFA [27]. The qualitative version of the FFA is identical to the regular FFA except that it only tests the undiluted sample and therefore does not allow end-point titration. For a small subset of samples (S2 File) the FFA result was not reliable because of fungal contamination of the well. These undetermined samples were re-tested by two-step reverse transcription polymerase chain reaction (RT-PCR) detection of DENV RNA as previously described [27].

**Virus sequencing.** The full-length genomes of two of the DENV isolates (Thai2010a and Thai2010b) were already sequenced from a previous study [41] and are available in GenBank under accession numbers HG316481 and HG316482, respectively. The full-length genome sequences of the six other DENV isolates of the study were obtained by high-throughput sequencing. Total RNA was extracted from cell culture supernatant using QIAamp Viral RNA Mini Kit (Qiagen), DNAse treated (Turbo DNAse, Life Technologies) and purified with magnetic beads (Beckman Coulter). Purified RNA was then reverse transcribed using Transcription High Fidelity cDNA Synthesis Kit and random hexamers (Roche Applied Science, Penzberg, Germany). Second DNA strand was synthetized in a single reaction with *E. coli* DNA ligase (NEB), *E. coli* DNA polymerase I (NEB), *E. coli* RNAse H (NEB) in second strand synthesis buffer (NEB). Resulting dsDNA was purified with magnetic beads (Agencourt AMPure XP, Beckman Coulter) and its concentration was measured by fluorometric quantification (Quant-iT PicoGreen dsDNA, Invitrogen). Sequencing libraries were prepared using Nextera XT DNA Library Preparation Kit (Illumina), multiplexed and sequenced in single-end on an Illumina NextSeq 500 platform using a mid-output 150-cycle v2 kit (Illumina).

After demultiplexing, trimmomatic v0.33 [42] was used to discard reads shorter than 32 nucleotides, filter out Illumina adaptor sequences, remove leading and trailing low quality bases and trim reads when the average quality per base dropped below 15 on a 4-base-wide sliding window. Full DENV genomes were reconstructed by *de novo* assembly using Ray v2.0.0 [43]. The largest contig was extended in 3' and 5' using the full genome sequence of the closest BLAST hit. This chimeric construct was used as the reference to re-map all the reads using Bowtie 2 v2.1.0 [44] and correct the 3' and 5' ends. The alignment file was converted, sorted and indexed using Samtools v0.1.19 [45]. Coverage and sequencing depth was assessed for each sample using bedtools v2.17.0 [46]. Single nucleotide variants and their frequency were called using LoFreq\* v2.1.1 [47] and used to correct the chimeric construct. Only nucleotides with >10X coverage were conserved for generating the consensus sequence. The full-length genomes of the six newly sequenced DENV isolates were deposited in GenBank under accessions numbers MG877552 (Laos2012), MG877553 (Haiti2012), MG877554 (Thai2012), MG877555 (Ncal2013), MG877556 (Thai2013), and MG877557 (Gabon2012).

**Phylogenetic analyses.** A background set of 51 full-length DENV genome sequences was obtained from GenBank to represent the different genotypes of DENV type 1 [48]. Genome sequences were aligned using ClustalW v2.0.12 [49] and analyzed with RAxML v8.2.10 [50] to generate the best-scoring maximum-likelihood (ML) tree of 10 runs with 100 thorough bootstrap replicates. The GTR + I + G nucleotide substitution model was chosen based on the lowest Bayesian Information Criterion (BIC) value using the JModelTest 2 software [51]. The same procedure was applied to the eight DENV isolates. Phylogenetic trees were visualized using FigTree v1.4.3 [52] and the ggtree R package [53]. Phylogenetic comparative analysis of continuous traits was performed with the phyloSignal function of the phylosignal R package [54] using the Abouheif's Cmean, Moran's I, Blomberg's K and Pagel's Lambda methods [55].

**Empirical data analysis.** A global likelihood function was used to model the probabilities of systemic mosquito infection over time for each virus isolate. Probabilities of systemic
infection at each time point post virus exposure were estimated with a 3-parameter logistic model (Eq 1) assuming a binomial distribution of systemic infection status.

$$f(t, (K, B, M)) = K/(1 + e^{-B(t-M)}) \qquad \text{Eq 1}$$

Based on this function, the probability of systemic infection at a given time point (t) is governed by three parameters. K is the saturation level and represents the maximum proportion of mosquitoes with a systemic infection. B is the slope factor and represents the maximum value of the slope during the exponential phase of the cumulative function, scaled by K. For easier biological interpretation, B was transformed into Δt, which corresponds to the time required to rise from 10% to 90% of the saturation level with the formula: Δt = ln(81)/B [56, 57]. M is the lag time and represents the time at which the absolute increase in cumulative proportion is maximal. For each isolate, the subplex R function [58] was used to provide the best estimates of the three parameters to maximize the global likelihood function (i.e., the sum of binomial probabilities at each time point post virus exposure). This method accounts for differences in sample size when estimating parameters values. The daily rate of new systemic infections was obtained by calculating the first derivative of the cumulative systemic infection function.

Comparison of systemic infection kinetics among DENV isolates was based on the Akaike Information Criterion (AIC). The global likelihood model was fitted to all 4,140 partitions of isolates. Isolates from one group were constrained to share the same parameter values. AIC values were calculated for all possible partitions of isolates with the same empirical data of systemic infection. The best isolate grouping was considered as the one with the lowest AIC score. An affinity matrix among isolates was derived from AIC values obtained from all grouping comparisons. Two isolates share a higher affinity if they are commonly grouped together in models with low AIC values. A hierarchical clustering analysis was performed on the affinity scores to reveal the extent of similarities of systemic infection kinetics among isolates.

The time-dependent effect of infectious dose (blood meal titer) on mosquito midgut infection and systemic infection was analyzed by logistic regression. Differences in mosquito infection over time were analyzed using hierarchical clustering with the Ward method based on the scaled proportion of infected mosquitoes at each time point post exposure. Pvclust [59] was used to calculate Approximately Unbiased (AU) p values for hierarchical clustering through multiscale bootstrap resampling using 1,000 bootstrap replicates. Genome-wide associations were tested by analysis of variance of the kinetics parameter values as a function of viral sequence polymorphisms (at the nucleotide and amino-acid levels). The null hypothesis was the lack of difference in parameter values between polymorphic variants at each locus. It was rejected when p values were <0.05 after Bonferroni correction for multiple testing. Geneious v.7.0.6 (http://www.geneious.com, [60]) was used to align and translate the viral genomic sequences. All statistical analyses were performed in the R statistical environment [61] and the R Markdown in provided as S3 File. Charts were generated with the ggplot2 package [62].

## Epidemiological modeling

**Model structure overview.** A stochastic agent-based model was developed in a closed human population to assess the influence of variation in the dynamics of systemic mosquito infection on dengue outbreaks. An agent-based model was chosen to achieve a more intuitive integration of the complex interplay of relatively simple rules within and between agents than would be possible with a mechanistic model. Starting with one infected human, this model encapsulates arbovirus transmission between susceptible mosquitoes and humans within a simulated environment. The model considers the influence of intra-host infection dynamics,
both in humans and mosquitoes, on virus transmission probability during human-mosquito contacts. The model was run for several combinations of the three parameters (K, B and M) that were used above to describe the kinetics of systemic mosquito infection. The main epidemiological outcome was the total number of secondary human infections. The entities of the model were individual humans and mosquitoes. The variable states of humans were uninfected/infected. The variable states of mosquitoes were alive/dead and uninfected/infected. The time step of the model was one day. The values assigned to the model variables are shown in S1 Table.

**Human and vector populations.** Human and mosquito population sizes were set as constant for each simulation. The model was populated with 10,000 humans and 30,000 mosquitoes. No birth, death or migration was allowed for humans. At each time point, daily mosquito deaths were compensated by the birth of uninfected mosquitoes. The model assumed homogeneous mixing of humans and mosquitoes, allowing any mosquito to bite any human with equal probability. At each time point in the simulation, a subset of the mosquito population was randomly chosen based on the daily feeding rate (S1 Table) to bite randomly chosen humans. The maximal number of daily bites of an individual mosquito was limited to one, whereas several mosquitoes could bite an individual human in one day.

**Transmission event.** The model only considered arbovirus transmission (i) from an infected human to uninfected vectors and (ii) from an infected vector to susceptible humans. The model ignored vertical transmission from an infected female to its offspring because it is rare both in humans [63] and in mosquitoes [64] and can be considered as epidemiologically insignificant. The model assumed that humans acquired protective immunity from a previous infection, whereas infection clearance was not allowed in mosquitoes. No superinfection events were allowed in mosquitoes and humans that would change infection dynamics.

**Human-to-mosquito transmission.** Dynamics of human-to-mosquito transmission probability was modeled using a 5-parameter equation (Eq 2) where the transmission probability *pvert* is linked to time since infection (t), the maximum probability of transmission (*K*<sub>human</sub>), time when 50% of *K*<sub>human</sub> value is attained during ascending probability (*M*γ), slope of the ascending curve (*B*γ), time when 50% of *K*<sub>human</sub> value is attained during descending probability (*M*δ), and slope of the descending curve (*B*δ).

$$pvert(t) = K_{\text{human}}/(1 + e^{(-B\gamma \times (t-M\gamma))}) - K_{\text{human}}/(1 + e^{(-B\delta \times (t-M\delta))}) \quad \text{Eq 2}$$

For each individual human, values for these five parameters were drawn from a probability distribution, of which specific parameters are provided in S1 Table. *M*γ, *B*γ, and *B*δ were drawn from a normal probability distribution, *δ* was drawn from a truncated normal probability distribution (to avoid values <*M*γ), and *K*<sub>human</sub> was drawn from a beta distribution. Distributions of each parameter were set based on previously published data on DENV infection dynamics in humans (S1 Table). By reanalyzing classic studies in which DENV-infected human volunteers were experimentally exposed to mosquitoes [19, 25], Nishiura and Halstead inferred probabilities of successful transmission from humans to mosquitoes according to time of disease onset in humans [65]. The link between time since infection and transmission probability was obtained by transforming time of disease onset in humans into time post human infection based on the work of Chan and Johansson [66]. Using a meta-analysis, Chan and Johansson estimated that the average time between infection and onset of the symptoms (the intrinsic incubation period) was about 5.9 days for dengue. Eq 1 was fitted to the transformed data to estimate parameter values presented in S1 Table.

**Mosquito-to-human transmission.** Dynamics of mosquito-to-human transmission probability was modeled using the 3-parameter equation (Eq 1) described above. For each
individual mosquito, values for the three parameters were drawn from a probability distribu-
tion (S1 Table). M was drawn from a normal probability distribution, B was drawn from a
truncated normal distribution (to avoid values <0) and K was drawn from a beta distribution.
Alternative modeling strategies are described in S1 File.

**Implementation.** The model was implemented in R version 3.4.0 using the packages fore-
ach [67], doParallel [68], plyr [69], reshape2 [70], and data.table [71]. A total of 1,140 different
parameter combinations for *K*, *B*, and *M* were run. Two sets of simulations were run. The first
set aimed to explore the parameter space of *K*, *B*, and *M*. An arbitrary range of mean values
was generated for each parameter by increasing from a realistic minimum value to a realistic
maximum value in a stepwise fashion (the increment size is provided in S1 Table). A total of
1,140 different parameter combinations for *K*, *B*, and *M* were run, and each combination was
run for 100 independent replicate simulations, yielding a total number of 114,000 simulations.
The second set of simulations aimed to explore epidemiological outcomes for the eight empiri-
cally characterized DENV isolates. In these simulations, mean values of *K*, *B*, and *M* were set
to their respective experimentally determined estimates, whereas standard deviation was the
same as in the first set of simulations (S1 Table). Epidemics caused by each DENV isolate were
examined in 100 independent replicate simulations. Each simulation was initiated with one
infected human at *t* = 0, then each simulation was run for 400 days. The code of simulations is
available at https://github.com/slequime/ArboEpiSim.

**In silico data analysis.** For each replicate of a given parameter combination, the total number of humans infected during the outbreak was recorded. Two outputs derived from the simulation data were examined: (i) the proportion of failures to initiate an outbreak and (ii) the mean number of secondary human infections across replicates. Failure to initiate an outbreak was defined as the lack of secondary infections over the length of the simulation following initial introduction of a single infected human. Failures to initiate an outbreak were filtered out before calculating the mean number of human infections. For each analysis, a stratified random splitting was performed on the data with the createDataPartition function from the caret R package [72] to obtain *training* and *test* subsets of the data with a balanced number of sample for each parameters combinations. The *training* set was composed of 70% of the data. We performed random forest regression with optimized parameter settings provided by the train function of the caret R package [72]. A total of 500 trees were used for both analyses with two and four variables or interaction between variables randomly tried at each split. The model included the proportion of failures to initiate an outbreak as the response variable and the three parameters describing systemic infection dynamics and their interaction as explanatory variables. The relative importance of parameters was determined according to their *IncNodePurity* value, which measures the quality of a split for each variable, averaged over all trees.

**Supporting information**

S1 Fig. Observed prevalence of mosquito midgut infection following oral exposure to DENV isolates. (A) Bar charts represent the proportion of mosquitoes with a DENV-infected midgut at several time points following oral exposure. Dashes represent the 95% confidence intervals of the proportions. Time-averaged proportions are shown on the right side of each chart. Midgut infection prevalence was significantly influenced by the time point (p = 6.6×10⁻⁹), and increased with uncontrolled variation in blood meal titers although not significantly so (p = 0.12). (B) The heat map represents midgut infection prevalence (using scaled values to standardize variation across time points) over time post DENV exposure for each isolate. Hierarchical clustering dendrograms are displayed on the left side of the heat map with
Approximately Unbiased (AU) values (in %) computed by multiscale bootstrap resampling. A green rectangle highlights isolates groupings with AU values larger than 95%. (EPS)

**S2 Fig. Pairwise correlations between the three parameters of systemic mosquito infection dynamics.** For each isolate, the cumulative change in the proportion of mosquitoes with a systemic DENV infection over time post exposure was used to fit a 3-parameter logistic model. *K* is the saturation level and represents the maximum proportion of mosquitoes with a systemic infection. *B* is the slope factor and represents the maximum slope of the cumulative function scaled by *K*. *M* is the lag time and represents the time at which the absolute increase in cumulative proportion is maximal. Blue lines represent linear regressions and the gray areas are their confidence intervals. Spearman's rank correlation coefficients and *p* values are indicated below each graph.
(EPS)

**S3 Fig. Statistical independence between parameter of systemic mosquito infection dynamics and phylogenetic relatedness among DENV isolates.** (A) Scaled and centered parameters (*K*: saturation level; *B*: slope factor; *M*: lag time) describing the kinetics of systemic mosquito infection are mapped along the phylogeny of the 8 DENV isolates. The tree was constructed with a maximum-likelihood method implemented in RAxML v.8.2.10 [50] using the full-genome sequence of the isolates. The scale bar represents the number of substitutions per site. (B) Statistical independence between parameter values and phylogenetic relatedness among isolates was assessed using the Abouheif's Cmean, Moran's I, Blomberg's K and Pagel's Lambda methods implemented in the phylosignal R package [54]. The table displays *p* values obtained with each method.
(EPS)

**S4 Fig. Analyses of simulations outputs.** (A) Density chart of the mean number of infected humans across simulation replicates for several combinations of *K* and *M* (with parameter *B* drawn from a uniform distribution). *K* and *M* were the two most influential parameters identified by random forest analysis performed on two independent simulation outputs: (B) the proportion of failures to initiate an outbreak and (C) the mean number of human infections across replicates. Overall, random forest analysis explained 88.9% and 96.6% of the variability in the proportion of failure to initiate an outbreak and in the mean number of human infections in successful outbreaks, respectively. In panels (B) and (C), upper charts show the correlation between observed and predicted values on the *test* data subset. Bottom charts show the relative importance of infection kinetics parameters, with higher *IncNode Purity* values meaning greater importance.
(EPS)

**S5 Fig. Identification of the best-fit model of systemic mosquito infection kinetics among all possible isolate groupings.** A logistic 3-parameter model was fitted to all permutations of isolates, with isolates from the same group being forced to share the same parameters. AIC values for each isolate grouping are represented as a function of the number of groups in the model. The inset shows AIC values for all isolate groupings. The main graph shows the isolate groupings in the lower range of AIC values. The lowest AIC value is indicated by a green circle and represents the isolate grouping that maximizes the probability function shown in Fig 3A.
(TIF)

**S1 Table. Central values used for each parameter in the agent-based model.** *Parameters drawn from a beta distribution. **Parameters drawn from a truncated normal distribution to
avoid non-existing values.

(DOCX)

**S1 File. Comparison of simulation outputs between EIP modeling strategies.** Simulations results are presented for the eight DENV isolates when EIP was modeled by a 3-parameter equation, a threshold effect with variable EIP, or a fixed EIP value (EIP$_{50}$ or EIP$_{10}$) for the entire mosquito population.

(PDF)

**S2 File. Raw empirical data from this study.** Each line represents an individual mosquito and the columns are the experimental conditions and phenotypes. From left to right, the columns indicate the DENV isolate, number of days post infectious blood meal, body infection status, head infection status, infection round (out of three sequential 15-min rounds of exposure), experiment (out of two separate experiments), blood meal titer (in log$_{10}$ FFU/ml), number of prior passages in C6/36 cells, country of origin, date of isolation, and DENV-1 genotype. The last column identifies the small number of samples that were tested by PCR instead of FFA.

(TXT)

**S3 File. R Markdown document of the study.** The file provides code lines used to analyze the experimental data and perform data visualization. The R code lines only refer to the procedure of systemic infection modeling and the AIC-based method to compare systemic infection across isolates. This document can be opened in Rstudio [73] with S2 File as the input file.

(RMD)

**S4 File. Simulated effect of mosquito infection kinetics on dengue outbreaks.** The interactive bubble plot represents the probability and size of dengue outbreaks across simulation replicates as a function of the three parameters (K, B, M) describing the kinetics of systemic mosquito infection. The size of the bubble is proportional to the percentage of successful outbreaks (>1 infected human). The color of the bubble represents the log$_{10}$-transformed mean number of infected humans during successful outbreaks. Hover text indicates the percentage of successful outbreaks, the mean number of infected humans during successful outbreaks, and the parameter values when the pointing device is moved over a bubble. This chart was made with Plotly for R v. 4.7.1 [74] and can be viewed in most web browsers.

(HTML)

**Acknowledgments**

We are grateful to Alongkot Ponlawat and Thanyalak Fansiri for the initial field collection of mosquitoes in Thailand, and to Catherine Lallemand for assistance with mosquito rearing. We thank the Lambrechts lab members for helpful discussions, and three anonymous reviewers for constructive criticism on earlier versions of the manuscript.

**Author Contributions**

**Conceptualization:** Albin Fontaine, Sebastian Lequime, Louis Lambrechts.

**Formal analysis:** Albin Fontaine, Sebastian Lequime, Robert Charles Reiner, Jr, Louis Lambrechts.

**Funding acquisition:** Louis Lambrechts.

**Investigation:** Albin Fontaine, Sebastian Lequime, Isabelle Moltini-Conclois, Davy Jiolle.
**Methodology:** Albin Fontaine, Sebastian Lequime, Robert Charles Reiner, Jr, Louis Lambrechts.

**Resources:** Isabelle Leparc-Goffart.

**Supervision:** Louis Lambrechts.

**Visualization:** Albin Fontaine, Sebastian Lequime.

**Writing – original draft:** Albin Fontaine, Louis Lambrechts.

**Writing – review & editing:** Albin Fontaine, Sebastian Lequime, Isabelle Moltini-Conclois, Davy Jiolle, Robert Charles Reiner, Jr, Louis Lambrechts.

## References

1. Wilder-Smith A, Gubler DJ, Weaver SC, Monath TP, Heymann DL, Scott TW. Epidemic arboviral diseases: priorities for research and public health. Lancet Infect Dis. 2017; 17(3):e101-e6. https://doi.org/10.1016/S1473-3099(16)30518-7 PMID: 28011234.

2. Tabachnick WJ. Ecological effects on arbovirus-mosquito cycles of transmission. Curr Opin Virol. 2016; 21:124-31. https://doi.org/10.1016/j.coviro.2016.09.008 PMID: 27693828.

3. Kramer LD, Ciota AT. Dissecting vectorial capacity for mosquito-borne viruses. Curr Opin Virol. 2015; 15:112-8. https://doi.org/10.1016/j.coviro.2015.10.003 PMID: 26569343.

4. Perkins TA, Scott TW, Le Menach A, Smith DL. Heterogeneity, mixing, and the spatial scales of mosquito-borne pathogen transmission. PLoS Comput Biol. 2013; 9(12):e1003327. https://doi.org/10.1371/journal.pcbi.1003327 PMID: 24348223.

5. Streicker DG, Fenton A, Pedersen AB. Differential sources of host species heterogeneity influence the transmission and control of multihost parasites. Ecol Lett. 2013; 16(8):975-84. https://doi.org/10.1111/ele.12122 PMID: 23714379.

6. Vazquez-Prokopec GM, Perkins TA, Waller LA, Lloyd AL, Reiner RC Jr., Scott TW, et al. Coupled Heterogeneities and Their Impact on Parasite Transmission and Control. Trends Parasitol. 2016; 32 (5):356-67. https://doi.org/10.1016/j.pt.2016.01.001 PMID: 26850821.

7. Woolhouse ME, Dye C, Etard JF, Smith T, Charlwood JD, Garnett GP, et al. Heterogeneities in the transmission of infectious agents: implications for the design of control programs. Proc Natl Acad Sci U S A. 1997; 94(1):338-42. PMID: 8990210.

8. Smith DL, Battle KE, Hay SI, Barker CM, Scott TW, McKenzie FE. Ross, macdonald, and a theory for the dynamics and control of mosquito-transmitted pathogens. PLoS Pathog. 2012; 8(4):e1002588. https://doi.org/10.1371/journal.ppat.1002588 PMID: 22496640.

9. Reiner RC Jr., Perkins TA, Barker CM, Niu T, Chaves LF, Ellis AM, et al. A systematic review of mathematical models of mosquito-borne pathogen transmission: 1970-2010. J R Soc Interface. 2013; 10 (81):20120921. https://doi.org/10.1098/rsif.2012.0921 PMID: 23407571.

10. Lambrechts L, Paaijmans KP, Fansiri T, Carrington LB, Kramer LD, Thomas MB, et al. Impact of daily temperature fluctuations on dengue virus transmission by *Aedes aegypti*. Proc Natl Acad Sci U S A. 2011; 108(18):7460-5. https://doi.org/10.1073/pnas.1101377108 PMID: 21502510.

11. Liu-Helmersson J, Quam M, Wilder-Smith A, Stenlund H, Ebi K, Massad E, et al. Climate Change and *Aedes* Vectors: 21st Century Projections for Dengue Transmission in Europe. EBioMedicine. 2016; 7:267-77. https://doi.org/10.1016/j.ebiom.2016.03.046 PMID: 27322480.

12. Tjaden NB, Thomas SM, Fischer D, Beierkuhnlein C. Extrinsic Incubation Period of Dengue: Knowledge, Backlog, and Applications of Temperature Dependence. PLoS Negl Trop Dis. 2013; 7(6):e2207. https://doi.org/10.1371/journal.pntd.0002207 PMID: 23826399.

13. Bhatt S, Gething PW, Brady OJ, Messina JP, Farlow AW, Moyes CL, et al. The global distribution and burden of dengue. Nature. 2013; 496(7446):504-7. https://doi.org/10.1038/nature12060 PMID: 23563266.

14. Katzelnick LC, Fonville JM, Gromowski GD, Bustos Arriaga J, Green A, James SL, et al. Dengue viruses cluster antigenically but not as discrete serotypes. Science. 2015; 349(6254):1338-43. https://doi.org/10.1126/science.aac5017 PMID: 26383952.

15. Chen R, Vasilakis N. Dengue—quo tu et quo vadis? Viruses. 2011; 3(9):1562-608. https://doi.org/10.3390/v3091562 PMID: 21994796.

16. Lambrechts L, Scott TW, Gubler DJ. Consequences of the expanding global distribution of *Aedes albo pictus* for dengue virus transmission. PLoS Negl Trop Dis. 2010; 4(5):e646. https://doi.org/10.1371/journal.pntd.0000646 PMID: 20520794.
17. Lambrechts L, Fansiri T, Pongsiri A, Thaisomboonsuk B, Klungthong C, Richardson JH, et al. Dengue-1 virus clade replacement in Thailand associated with enhanced mosquito transmission. J Virol. 2012; 86(3):1853-61. https://doi.org/10.1128/JVI.06458-11 PMID: 22130539.

18. Christofferson RC, Mores CN. Estimating the magnitude and direction of altered arbovirus transmission due to viral phenotype. PLoS ONE. 2011; 6(1):e16298. https://doi.org/10.1371/journal.pone.0016298 PMID: 21298018.

19. Siler JF, Hall MW, Hitchens AP. Dengue: its history, epidemiology, mechanism of transmission, etiology, clinical manifestations, immunity, and prevention. Philippine J Sci. 1926; 29:1-302.

20. Watts DM, Burke DS, Harrison BA, Whitmire RE, Nisalak A. Effect of temperature on the vector efficiency of *Aedes aegypti* for dengue 2 virus. Am J Trop Med Hyg. 1987; 36(1):143-52. PMID: 3812879.

21. Carrington LB, Tran BCN, Le NTH, Luong TTH, Nguyen TT, Nguyen PT, et al. Field- and clinically derived estimates of *Wolbachia*-mediated blocking of dengue virus transmission potential in *Aedes aegypti* mosquitoes. Proc Natl Acad Sci U S A. 2018; 115(2):361-6. https://doi.org/10.1073/pnas.1715788115 PMID: 29279375.

22. Ferguson NM, Kien DT, Clapham H, Aguas R, Trung VT, Chau TN, et al. Modeling the impact on virus transmission of *Wolbachia*-mediated blocking of dengue virus infection of *Aedes aegypti*. Sci Transl Med. 2015; 7(279):279ra37. https://doi.org/10.1126/scitranslmed.3010370 PMID: 25787763.

23. Anderson JR, Rico-Hesse R. *Aedes aegypti* vectorial capacity is determined by the infecting genotype of dengue virus. Am J Trop Med Hyg. 2006; 75(5):886-92. PMID: 17123982.

24. Ohm JR, Baldini F, Barreaux P, Lefevre T, Lynch PA, Suh E, et al. Rethinking the extrinsic incubation period of malaria parasites. Parasites Vectors. 2018; 11(1):178. https://doi.org/10.1186/s13071-018-2761-4 PMID: 29530073.

25. Simmons JS, St. John JH, Reynolds FHK. Experimental studies of dengue. Philippine J Sci. 1931; 44:1-247.

26. Moudy RM, Meola MA, Morin LL, Ebel GD, Kramer LD. A newly emergent genotype of West Nile virus is transmitted earlier and more efficiently by *Culex* mosquitoes. Am J Trop Med Hyg. 2007; 77(2):365-70. 77/2/365 [pii]. PMID: 17690414.

27. Fontaine A, Jiolle D, Moltini-Conclois I, Lequime S, Lambrechts L. Excretion of dengue virus RNA by *Aedes aegypti* allows non-destructive monitoring of viral dissemination in individual mosquitoes. Sci Rep. 2016; 6:24885. https://doi.org/10.1038/srep24885 PMID: 27117953.

28. Paaijmans KP, Blanford S, Chan BH, Thomas MB. Warmer temperatures reduce the vectorial capacity of malaria mosquitoes. Biol Lett. 2012; 8(3):465-8. https://doi.org/10.1098/rsbl.2011.1075 PMID: 22188673.

29. Carrington LB, Armijos MV, Lambrechts L, Scott TW. Fluctuations at a low mean temperature accelerated dengue virus transmission by *Aedes aegypti*. PLoS Negl Trop Dis. 2013; 7(4):e2190. https://doi.org/10.1371/journal.pntd.0002190 PMID: 23638208.

30. Ye YH, Chenoweth SF, Carrasco AM, Allen SL, Frentiu FD, van den Hurk AF, et al. Evolutionary potential of the extrinsic incubation period of dengue virus in *Aedes aegypti*. Evolution. 2016; 70(11):2459-69. https://doi.org/10.1111/evo.13039 PMID: 27530960.

31. Mammen MP, Pimgate C, Koenraadt CJ, Rothman AL, Aldstadt J, Nisalak A, et al. Spatial and temporal clustering of dengue virus transmission in Thai villages. PLoS Med. 2008; 5(11):e205. https://doi.org/10.1371/journal.pmed.0050205 PMID: 18986209.

32. Lambrechts L, Scott TW. Mode of transmission and the evolution of arbovirus virulence in mosquito vectors. Proc Biol Sci. 2009; 276(1660):1369-78. https://doi.org/10.1098/rspb.2008.1709 PMID: 19141420.

33. Dubrulle M, Mousson L, Moutailler S, Vazeille M, Failloux AB. Chikungunya virus and *Aedes* mosquitoes: saliva is infectious as soon as two days after oral infection. PLoS ONE. 2009; 4(6):e5895. https://doi.org/10.1371/journal.pone.0005895 PMID: 19521520.

34. Scott TW, Burrage TG. Rapid infection of salivary glands in *Culiseta melanura* with eastern equine encephalitis virus: an electron microscopic study. Am J Trop Med Hyg. 1984; 33(5):961-4. PMID: 6148898.

35. Ng WC, Soto-Acosta R, Bradrick SS, Garcia-Blanco MA, Ooi EE. The 5' and 3' Untranslated Regions of the Flaviviral Genome. Viruses. 2017; 9(6). https://doi.org/10.3390/v9060137 PMID: 28587300.

36. Salazar MI, Richardson JH, Sanchez-Vargas I, Olson KE, Beaty BJ. Dengue virus type 2: replication and tropisms in orally infected *Aedes aegypti* mosquitoes. BMC Microbiol. 2007; 7:9. https://doi.org/10.1186/1471-2180-7-9 PMID: 17263893.

37. Black WC, Bennett KE, Gorrochotegui-Escalante N, Barillas-Mury CV, Fernandez-Salas I, de Lourdes Munoz M, et al. Flavivirus susceptibility in *Aedes aegypti*. Arch Med Res. 2002; 33(4):379-88. PMID: 12234528.
38. Duong V, Lambrechts L, Paul RE, Ly S, Lay RS, Long KC, et al. Asymptomatic humans transmit dengue virus to mosquitoes. Proc Natl Acad Sci U S A. 2015; 112(47):14688-93. https://doi.org/10.1073/pnas.1508114112 PMID: 26553981.

39. Nguyen N, Thi Hue Kien D, Tuan T, Quyen N, Tran C, Vo Thi L, et al. Host and viral features of human dengue cases shape the population of infected and infectious *Aedes aegypti* mosquitoes. Proc Natl Acad Sci U S A. 2013; 110(22):9072-7. https://doi.org/10.1073/pnas.1303395110 PMID: 23674683

40. Siettos CI, Russo L. Mathematical modeling of infectious disease dynamics. Virulence. 2013; 4(4):295-306. https://doi.org/10.4161/viru.24041 PMID: 23552814.

41. Fansiri T, Fontaine A, Diancourt L, Caro V, Thaisomboonsuk B, Richardson JH, et al. Genetic Mapping of Specific Interactions between *Aedes aegypti* Mosquitoes and Dengue Viruses. PLoS Genet. 2013; 9(8):e1003621. https://doi.org/10.1371/journal.pgen.1003621 PMID: 23935524.

42. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014; 30(15):2114-20. https://doi.org/10.1093/bioinformatics/btu170 PMID: 24695404.

43. Boisvert S, Raymond F, Godzaridis E, Laviolette F, Corbeil J. Ray Meta: scalable de novo metagenome assembly and profiling. Genome Biol. 2012; 13(12):R122. https://doi.org/10.1186/gb-2012-13-12-r122 PMID: 23259615.

44. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012; 9(4):357-9. https://doi.org/10.1038/nmeth.1923 PMID: 22388286.

45. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009; 25(16):2078-9. https://doi.org/10.1093/bioinformatics/btp352 PMID: 19505943.

46. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010; 26(6):841-2. https://doi.org/10.1093/bioinformatics/btq033 PMID: 20110278.

47. Wilm A, Aw PP, Bertrand D, Yeo GH, Ong SH, Wong CH, et al. LoFreq: a sequence-quality aware, ultra-sensitive variant caller for uncovering cell-population heterogeneity from high-throughput sequencing datasets. Nucleic Acids Res. 2012; 40(22):11189-201. https://doi.org/10.1093/nar/gks918 PMID: 23066108.

48. Pyke AT, Moore PR, Taylor CT, Hall-Mendelin S, Cameron JN, Hewitson GR, et al. Highly divergent dengue virus type 1 genotype sets a new distance record. Sci Rep. 2016; 6:22356. https://doi.org/10.1038/srep22356 PMID: 26924208.

49. Thompson JD, Higgins DG, Gibson TJ. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 1994; 22(22):4673-80. PMID: 7984417.

50. Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014; 30(9):1312-3. https://doi.org/10.1093/bioinformatics/btu033 PMID: 24451623.

51. Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012; 9(8):772. https://doi.org/10.1038/nmeth.2109 PMID: 22847109.

52. Rambaut A. FigTree v1.4.3, Available from http://tree.bio.ed.ac.uk/software/figtree/. 1.4.3 ed2016.

53. Yu G, Smith DK, Zhu H, Guan Y, Lam TTY. ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol. 2017; 8(1):28-36.

54. Keck F, Rimet F, Bouchez A, Franc A. phylosignal: an R package to measure, test, and explore the phylogenetic signal. Ecol Evol. 2016; 6(9):2774-80. https://doi.org/10.1002/ece3.2051 PMID: 27066252.

55. Revell LJ, Harmon LJ, Collar DC. Phylogenetic signal, evolutionary process, and rate. Syst Biol. 2008; 57(4):591-601. https://doi.org/10.1080/10635180802302427 PMID: 18709597.

56. Meyer PS, Ausubel JH. Carrying capacity: a model with logistically varying limits. Technol Forecast Soc Change. 1999; 61(3):209-14. https://doi.org/10.1016/S0040-1625(99)00022-0.

57. Meyer PS, Yung JW, Ausubel JH. A primer on logistic growth and substitution: the mathematics of the loglet lab software. Technol Forecast Soc Change. 1999; 61(3):247-71. https://doi.org/10.1016/S0040-1625(99)00021-9.

58. King AA. subplex: unconstrained optimization using the subplex algorithm. R package version 3.4.0; 2015.

59. Suzuki R, Shimodaira H. Pvclust: an R package for assessing the uncertainty in hierarchical clustering. R package version 3.4.0; 2006.

60. Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012; 28(12):1647-9. https://doi.org/10.1093/bioinformatics/bts199 PMID: 22543367.
61. R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2014.

62. Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer; 2009.

63. Chye JK, Lim CT, Ng KB, Lim JMH, George R, Lam SK. Vertical transmission of dengue. Clin Infect Dis. 1997; 25(6):1374-7. doi: 10.1086/516126. PMID: 9431381

64. Lequime S, Paul RE, Lambrechts L. Determinants of Arbovirus Vertical Transmission in Mosquitoes. PLoS Pathog. 2016; 12(5):e1005548. https://doi.org/10.1371/journal.ppat.1005548 PMID: 27171170.

65. Nishiura H, Halstead SB. Natural history of dengue virus (DENV)-1 and DENV-4 infections: reanalysis of classic studies. J Infect Dis. 2007; 195(7):1007-13. https://doi.org/10.1086/511825 PMID: 17330791.

66. Chan M, Johansson MA. The incubation periods of Dengue viruses. PLoS ONE. 2012; 7(11):e50972. https://doi.org/10.1371/journal.pone.0050972 PMID: 23226436.

67. Analytics R, Weston S. foreach: provides foreach looping construct for R. R package 3.4.0; 2015.

68. Analytics R, Weston S. doParallel: foreach parallel adaptor for the parallel package. R package version 3.4.0; 2014.

69. Wickham H. The Split-Apply-Combine Strategy for Data Analysis. J Stat Softw. 2011; 40(1):1-29. WOS:000289228100001.

70. Wickham H. Reshaping data with the reshape package. J Stat Softw. 2007; 21(12):1-20. WOS:000252429400001.

71. Dowle M, Short T, Lianoglou S. data.table: extension of data.frame for fast indexing, fast ordered joins, fast assignment, fast grouping and list columns. R package version 3.4.0; 2013.

72. Kuhn M. Caret: classification and regression training. R package version 3.4.0; 2015.

73. RStudio Team. RStudio: Integrated Development for R. Boston, MA: RStudio, Inc.; 2015.

74. Sievert C. plotly: create interactive web graphics via Plotly's javascript graphing library. R package version 3.4.0; 2016.